你已经派生过 gprMax
镜像自地址
https://gitee.com/sunhf/gprMax.git
已同步 2025-08-07 15:10:13 +08:00
Formatting cleanups
这个提交包含在：
Craig Warren
@@ -963,6 +963,8 @@ class DiscretePlaneWave(UserObjectMulti):
|
||||
p1: tuple required for the upper right position (x, y, z) of the total
|
||||
field, scattered field (TFSF) box.
|
||||
waveform_id: string required for identifier of waveform used with source.
|
||||
material_id: string optional of material identifier to use as the
|
||||
background material in the TFSF box.
|
||||
start: float optional to delay start time (secs) of source.
|
||||
stop: float optional to time (secs) to remove source.
|
||||
"""
|
||||
|
||||
@@ -25,8 +25,14 @@ import gprMax.config as config
|
||||
from .fields_outputs import Ix, Iy, Iz
|
||||
from .utilities.utilities import round_value
|
||||
|
||||
import scipy.constants as constants
|
||||
from .cython.planeWaveModules import getIntegerForAngles, getProjections, calculate1DWaveformValues, updatePlaneWave, getSource
|
||||
from .cython.planeWaveModules import (
|
||||
getIntegerForAngles,
|
||||
getProjections,
|
||||
calculate1DWaveformValues,
|
||||
updatePlaneWave,
|
||||
getSource,
|
||||
)
|
||||
|
||||
|
||||
class Source:
|
||||
"""Super-class which describes a generic source."""
|
||||
@@ -647,362 +653,533 @@ class TransmissionLine(Source):
|
||||
|
||||
|
||||
class DiscretePlaneWave(Source):
|
||||
"""Implements the discrete plane wave (DPW) formulation as described in
|
||||
Tan, T.; Potter, M. (2010). FDTD Discrete Planewave (FDTD-DPW)
|
||||
Formulation for a Perfectly Matched Source in TFSF Simulations., 58(8),
|
||||
0–2648. doi:10.1109/tap.2010.2050446
|
||||
"""
|
||||
Class to implement the discrete plane wave (DPW) formulation as described in
|
||||
Tan, T.; Potter, M. (2010).
|
||||
FDTD Discrete Planewave (FDTD-DPW) Formulation for a Perfectly Matched Source in TFSF Simulations. ,
|
||||
58(8), 0–2648. doi:10.1109/tap.2010.2050446
|
||||
__________________________
|
||||
|
||||
Instance variables:
|
||||
--------------------------
|
||||
m, int array : stores the integer mappings, m_x, m_y, m_z which determine the rational angles
|
||||
last element stores max(m_x, m_y, m_z)
|
||||
directions, int array : stores the directions of propagation of the DPW
|
||||
dimensions, int : stores the number of dimensions in which the simulation is run (2D or 3D)
|
||||
time_dimension, int : stores the time length over which the simulation is run
|
||||
E_fileds, double array : stores the electric flieds associated with the 1D DPW
|
||||
H_fields, double array : stores the magnetic fields associated with the 1D DPW
|
||||
"""
|
||||
|
||||
|
||||
def __init__(self, G):
|
||||
"""
|
||||
Defines the instance variables of class DiscretePlaneWave()
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
time_dimension, int : local variable to store the time length over which the simulation is run
|
||||
dimensions, int : local variable to store the number of dimensions in which the simulation is run
|
||||
def __init__(self, G):
|
||||
"""
|
||||
Args:
|
||||
m: int array stores the integer mappings, m_x, m_y, m_z which
|
||||
determine the rational angles last element stores
|
||||
max(m_x, m_y, m_z).
|
||||
directions: int array stores the directions of propagation of DPW.
|
||||
dimensions: int stores the number of dimensions in which the
|
||||
simulation is run (2D or 3D).
|
||||
time_dimension: int stores the time length over which the simulation
|
||||
is run.
|
||||
E_fields: double array stores the electric flieds associated with
|
||||
1D DPW.
|
||||
H_fields: double array stores the magnetic fields associated with
|
||||
1D DPW.
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
|
||||
super().__init__()
|
||||
self.m = np.zeros(3+1, dtype=np.int32) #+1 to store the max(m_x, m_y, m_z)
|
||||
self.m = np.zeros(3 + 1, dtype=np.int32) # +1 to store the max(m_x, m_y, m_z)
|
||||
self.directions = np.zeros(3, dtype=np.int32)
|
||||
self.length = 0
|
||||
self.projections = np.zeros(3, dtype=config.sim_config.dtypes["float_or_double"])
|
||||
self.projections = np.zeros(
|
||||
3, dtype=config.sim_config.dtypes["float_or_double"]
|
||||
)
|
||||
self.corners = None
|
||||
self.materialID = 1
|
||||
self.ds = 0
|
||||
|
||||
def initializeDiscretePlaneWave(self, psi, phi, Delta_phi, theta, Delta_theta, G):
|
||||
"""
|
||||
Method to create a DPW, assign memory to the grids and get field values at different time and space indices
|
||||
__________________________
|
||||
self.ds = 0
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
psi, float : polarization angle of the incident plane wave (in radians)
|
||||
phi, float : azimuthal angle of the incident plane wave (in radians)
|
||||
Delta_phi, float : permissible error in the rational angle approximation to phi (in radians)
|
||||
theta, float : polar angle of the incident plane wave (in radians)
|
||||
Delta_theta, float : permissible error in the rational angle approximation to theta (in radians)
|
||||
number, int : number of cells in the 3D FDTD simulation
|
||||
dx, double : separation between adjacent cells in the x direction
|
||||
dy, double : separation between adjacent cells in the y direction
|
||||
dz, double : separation between adjacent cells in the z direction
|
||||
dt, double : time step for the FDTD simulation
|
||||
__________________________
|
||||
def initializeDiscretePlaneWave(self, psi, phi, Delta_phi, theta, Delta_theta, G):
|
||||
"""Creates a DPW, assigns memory to the grids, and gets field values
|
||||
at different time and space indices.
|
||||
|
||||
Args:
|
||||
psi: float for polarization angle of the incident plane
|
||||
wave.
|
||||
phi: float for azimuthal angle (radians) of the incident plane wave.
|
||||
Delta_phi: float for permissible error in the rational angle
|
||||
(radians) approximation to phi.
|
||||
theta: float for polar angle (radians) of the incident plane wave.
|
||||
Delta_theta: float for permissible error in the rational angle
|
||||
(radians) approximation to theta.
|
||||
number: int for number of cells in the 3D FDTD simulation.
|
||||
dx: double for separation between adjacent cells in the x direction.
|
||||
dy: double for separation between adjacent cells in the y direction.
|
||||
dz: double for separation between adjacent cells in the z direction.
|
||||
dt: double for time step for the FDTD simulation.
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
E_fields, double array : the electric field for the DPW as it evolves over time and space indices
|
||||
H_fields, double array : the magnetic field for the DPW as it evolves over time and space indices
|
||||
C, double array : stores the coefficients of the fields for the update equation of the electric fields
|
||||
D, double array : stores the coefficients of the fields for the update equation of the magnetic fields
|
||||
|
||||
E_fields: double array for electric field for the DPW as it evolves
|
||||
over time and space indices.
|
||||
H_fields: double array for magnetic field for the DPW as it evolves
|
||||
over time and space indices.
|
||||
C: double array stores coefficients of the fields for the update
|
||||
equation of the electric fields.
|
||||
D: double array stores coefficients of the fields for the update
|
||||
equation of the magnetic fields.
|
||||
"""
|
||||
self.directions, self.m[:3] = getIntegerForAngles(phi, Delta_phi, theta, Delta_theta,
|
||||
np.array([G.dx, G.dy, G.dz])) #get the integers for the nearest rational angle
|
||||
self.m[3] = max(self.m) #store max(m_x, m_y, m_z) in the last element of the array
|
||||
self.length = int( 2 * max(self.m[:-1]) * G.iterations) #set an appropriate length fo the one dimensional arrays
|
||||
self.E_fields = np.zeros((3, self.length), order='C', dtype=config.sim_config.dtypes["float_or_double"])
|
||||
self.H_fields = np.zeros((3, self.length), order='C', dtype=config.sim_config.dtypes["float_or_double"])
|
||||
#the 1D grid has no ABC to terminate it, sufficiently long array prevents reflections from the back
|
||||
|
||||
self.directions, self.m[:3] = getIntegerForAngles(
|
||||
phi, Delta_phi, theta, Delta_theta, np.array([G.dx, G.dy, G.dz])
|
||||
) # get the integers for the nearest rational angle
|
||||
self.m[3] = max(
|
||||
self.m
|
||||
) # store max(m_x, m_y, m_z) in the last element of the array
|
||||
self.length = int(
|
||||
2 * max(self.m[:-1]) * G.iterations
|
||||
) # set an appropriate length fo the one dimensional arrays
|
||||
self.E_fields = np.zeros(
|
||||
(3, self.length),
|
||||
order="C",
|
||||
dtype=config.sim_config.dtypes["float_or_double"],
|
||||
)
|
||||
self.H_fields = np.zeros(
|
||||
(3, self.length),
|
||||
order="C",
|
||||
dtype=config.sim_config.dtypes["float_or_double"],
|
||||
)
|
||||
# the 1D grid has no ABC to terminate it, sufficiently long array prevents reflections from the back
|
||||
# Projections for field components
|
||||
projections_h, P = getProjections(psi*180/np.pi, self.m) #get the projection vertors for different fields
|
||||
self.projections = projections_h / np.sqrt(config.m0/config.e0) #scale the projection vector for the mangetic field
|
||||
|
||||
if self.m[0] == 0: #calculate dr that is needed for sourcing the 1D array
|
||||
projections_h, P = getProjections(
|
||||
psi * 180 / np.pi, self.m
|
||||
) # get the projection vertors for different fields
|
||||
self.projections = projections_h / np.sqrt(
|
||||
config.m0 / config.e0
|
||||
) # scale the projection vector for the mangetic field
|
||||
|
||||
if self.m[0] == 0: # calculate dr that is needed for sourcing the 1D array
|
||||
if self.m[1] == 0:
|
||||
if self.m[2] == 0:
|
||||
raise ValueError("not all m_i values can be zero")
|
||||
else:
|
||||
self.ds = P[2]*G.dz/self.m[2]
|
||||
self.ds = P[2] * G.dz / self.m[2]
|
||||
else:
|
||||
self.ds = P[1]*G.dy/self.m[1]
|
||||
self.ds = P[1] * G.dy / self.m[1]
|
||||
else:
|
||||
self.ds = P[0]*G.dx/self.m[0]
|
||||
self.ds = P[0] * G.dx / self.m[0]
|
||||
|
||||
|
||||
def calculate_waveform_values(self, G, cythonize=True):
|
||||
"""
|
||||
Calculates all waveform values for source for duration of simulation.
|
||||
"""Calculates all waveform values for source for duration of simulation.
|
||||
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
Args:
|
||||
G: FDTDGrid class describing a grid in a model.
|
||||
"""
|
||||
|
||||
# Waveform values for sources that need to be calculated on whole timesteps
|
||||
self.waveformvalues_wholedt = np.zeros((G.iterations, 3, self.m[3]), dtype=config.sim_config.dtypes["float_or_double"])
|
||||
self.waveformvalues_wholedt = np.zeros(
|
||||
(G.iterations, 3, self.m[3]),
|
||||
dtype=config.sim_config.dtypes["float_or_double"],
|
||||
)
|
||||
|
||||
# Waveform values for sources that need to be calculated on half timesteps
|
||||
#self.waveformvalues_halfdt = np.zeros((G.iterations), dtype=config.sim_config.dtypes["float_or_double"])
|
||||
# self.waveformvalues_halfdt = np.zeros((G.iterations), dtype=config.sim_config.dtypes["float_or_double"])
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
if(cythonize == True):
|
||||
calculate1DWaveformValues(self.waveformvalues_wholedt, G.iterations, self.m, G.dt, self.ds, config.c,
|
||||
self.start, self.stop, waveform.freq, waveform.type.encode('UTF-8'))
|
||||
if cythonize:
|
||||
calculate1DWaveformValues(
|
||||
self.waveformvalues_wholedt,
|
||||
G.iterations,
|
||||
self.m,
|
||||
G.dt,
|
||||
self.ds,
|
||||
config.c,
|
||||
self.start,
|
||||
self.stop,
|
||||
waveform.freq,
|
||||
waveform.type.encode("UTF-8"),
|
||||
)
|
||||
else:
|
||||
for dimension in range(3):
|
||||
for iteration in range(G.iterations):
|
||||
for iteration in range(G.iterations):
|
||||
for r in range(self.m[3]):
|
||||
time = G.dt * iteration - (r+ (self.m[(dimension+1)%3]+self.m[(dimension+2)%3])*0.5) * self.ds/config.c
|
||||
time = (
|
||||
G.dt * iteration
|
||||
- (
|
||||
r
|
||||
+ (
|
||||
self.m[(dimension + 1) % 3]
|
||||
+ self.m[(dimension + 2) % 3]
|
||||
)
|
||||
* 0.5
|
||||
)
|
||||
* self.ds
|
||||
/ config.c
|
||||
)
|
||||
if time >= self.start and time <= self.stop:
|
||||
# Set the time of the waveform evaluation to account for any
|
||||
# delay in the start
|
||||
time -= self.start
|
||||
self.waveformvalues_wholedt[iteration, dimension, r] = waveform.calculate_value(time, G.dt)
|
||||
#self.waveformvalues_halfdt[iteration] = waveform.calculate_value(time + 0.5 * G.dt, G.dt)
|
||||
self.waveformvalues_wholedt[iteration, dimension, r] = (
|
||||
waveform.calculate_value(time, G.dt)
|
||||
)
|
||||
# self.waveformvalues_halfdt[iteration] = waveform.calculate_value(time + 0.5 * G.dt, G.dt)
|
||||
|
||||
|
||||
def update_plane_wave(self, nthreads, updatecoeffsE, updatecoeffsH, Ex, Ey, Ez, Hx, Hy, Hz, iteration, G, cythonize=True, precompute=True):
|
||||
if(cythonize == True):
|
||||
def update_plane_wave(
|
||||
self,
|
||||
nthreads,
|
||||
updatecoeffsE,
|
||||
updatecoeffsH,
|
||||
Ex,
|
||||
Ey,
|
||||
Ez,
|
||||
Hx,
|
||||
Hy,
|
||||
Hz,
|
||||
iteration,
|
||||
G,
|
||||
cythonize=True,
|
||||
precompute=True,
|
||||
):
|
||||
if cythonize:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
updatePlaneWave(self.length, nthreads, self.H_fields, self.E_fields,
|
||||
updatecoeffsE[self.material_ID, :], updatecoeffsH[self.material_ID, :],
|
||||
Ex, Ey, Ez, Hx, Hy, Hz, self.projections, self.waveformvalues_wholedt[iteration, :, :],
|
||||
self.m, self.corners, precompute, iteration, G.dt, self.ds, config.c, self.start, self.stop,
|
||||
waveform.freq, waveform.type.encode('UTF-8'))
|
||||
updatePlaneWave(
|
||||
self.length,
|
||||
nthreads,
|
||||
self.H_fields,
|
||||
self.E_fields,
|
||||
updatecoeffsE[self.material_ID, :],
|
||||
updatecoeffsH[self.material_ID, :],
|
||||
Ex,
|
||||
Ey,
|
||||
Ez,
|
||||
Hx,
|
||||
Hy,
|
||||
Hz,
|
||||
self.projections,
|
||||
self.waveformvalues_wholedt[iteration, :, :],
|
||||
self.m,
|
||||
self.corners,
|
||||
precompute,
|
||||
iteration,
|
||||
G.dt,
|
||||
self.ds,
|
||||
config.c,
|
||||
self.start,
|
||||
self.stop,
|
||||
waveform.freq,
|
||||
waveform.type.encode("UTF-8"),
|
||||
)
|
||||
else:
|
||||
self.update_magnetic_field_1D(G, precompute)
|
||||
self.apply_TFSF_conditions_magnetic(G)
|
||||
self.apply_TFSF_conditions_electric(G)
|
||||
self.update_electric_field_1D(G)
|
||||
|
||||
|
||||
def initialize_magnetic_fields_1D(self, G, precompute):
|
||||
if(precompute == True):
|
||||
if precompute:
|
||||
for dimension in range(3):
|
||||
for r in range(self.m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
|
||||
self.H_fields[dimension, r] = self.projections[dimension] * self.waveformvalues_wholedt[G.iteration, dimension, r]
|
||||
#self.getSource(self.real_time - (j+(self.m[(i+1)%3]+self.m[(i+2)%3])*0.5)*self.ds/config.c)#, self.waveformID, G.dt)
|
||||
for r in range(
|
||||
self.m[3]
|
||||
): # loop to assign the source values of magnetic field to the first few gridpoints
|
||||
self.H_fields[dimension, r] = (
|
||||
self.projections[dimension]
|
||||
* self.waveformvalues_wholedt[G.iteration, dimension, r]
|
||||
)
|
||||
# self.getSource(self.real_time - (j+(self.m[(i+1)%3]+self.m[(i+2)%3])*0.5)*self.ds/config.c)#, self.waveformID, G.dt)
|
||||
else:
|
||||
waveform = next(x for x in G.waveforms if x.ID == self.waveformID)
|
||||
for dimension in range(3):
|
||||
for r in range(self.m[3]): #loop to assign the source values of magnetic field to the first few gridpoints
|
||||
self.H_fields[dimension, r] = self.projections[dimension] * getSource(G.iteration*G.dt - (r+(self.m[(dimension+1)%3]+self.m[(dimension+2)%3])*0.5)*self.ds/config.c, waveform.freq, waveform.type.encode('UTF-8'), G.dt)
|
||||
|
||||
for r in range(
|
||||
self.m[3]
|
||||
): # loop to assign the source values of magnetic field to the first few gridpoints
|
||||
self.H_fields[dimension, r] = self.projections[
|
||||
dimension
|
||||
] * getSource(
|
||||
G.iteration * G.dt
|
||||
- (
|
||||
r
|
||||
+ (
|
||||
self.m[(dimension + 1) % 3]
|
||||
+ self.m[(dimension + 2) % 3]
|
||||
)
|
||||
* 0.5
|
||||
)
|
||||
* self.ds
|
||||
/ config.c,
|
||||
waveform.freq,
|
||||
waveform.type.encode("UTF-8"),
|
||||
G.dt,
|
||||
)
|
||||
|
||||
def update_magnetic_field_1D(self, G, precompute=True):
|
||||
"""
|
||||
Method to update the magnetic fields for the next time step using
|
||||
Equation 8 of DOI: 10.1109/LAWP.2009.2016851
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
n , int : stores the spatial length of the DPW array so that each length grid cell is updated when the method updateMagneticFields() is called
|
||||
H_coefficients, double array : stores the coefficients of the fields in the update equation for the magnetic field
|
||||
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
|
||||
E_fields, double array : stores the electric fields of the DPW till temporal index time
|
||||
time, int : time index storing the current axis number which would be updated for the H_fields
|
||||
__________________________
|
||||
"""Updates magnetic fields for the next time step using Equation 8 of
|
||||
DOI: 10.1109/LAWP.2009.2016851
|
||||
|
||||
Args:
|
||||
n: int stores spatial length of the DPW array so that each length
|
||||
grid cell is updated when updateMagneticFields() called.
|
||||
H_coefficients: double array stores coefficients of the fields in
|
||||
the update equation for the magnetic field.
|
||||
H_fields: double array stores magnetic fields of the DPW until
|
||||
temporal index time.
|
||||
E_fields: double array stores electric fields of the DPW until
|
||||
temporal index time.
|
||||
time: int time index storing current axis number which would be
|
||||
updated for the H_fields.
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
H_fields, double array : magnetic field array with the axis entry for the current time added
|
||||
|
||||
H_fields: double array for magnetic field with the axis entry for
|
||||
the current time added.
|
||||
"""
|
||||
|
||||
self.initialize_magnetic_fields_1D(G, precompute)
|
||||
|
||||
for i in range(3): #loop to update each component of the magnetic field
|
||||
materialH = G.ID[3+i, (self.corners[0]+self.corners[3])//2,
|
||||
(self.corners[1]+self.corners[4])//2,
|
||||
(self.corners[2]+self.corners[5])//2 ]
|
||||
for j in range(self.m[-1], self.length-self.m[-1]): #loop to update the magnetic field at each spatial index
|
||||
self.H_fields[i, j] = G.updatecoeffsH[materialH, 0] * self.H_fields[i, j] + G.updatecoeffsH[materialH, (i+2)%3+1] * (
|
||||
self.E_fields[(i+1)%3, j+self.m[(i+2)%3]] - self.E_fields[(i+1)%3, j]) - G.updatecoeffsH[materialH, (i+1)%3+1] * (
|
||||
self.E_fields[(i+2)%3, j+self.m[(i+1)%3]] - self.E_fields[(i+2)%3, j]) #equation 8 of Tan, Potter paper
|
||||
|
||||
|
||||
def update_electric_field_1D(self, G):
|
||||
"""
|
||||
Method to update the electric fields for the next time step using
|
||||
Equation 9 of DOI: 10.1109/LAWP.2009.2016851
|
||||
__________________________
|
||||
|
||||
Input parameters:
|
||||
--------------------------
|
||||
n , int : stores the spatial length of the DPW array so that each length grid cell is updated when the method updateMagneticFields() is called
|
||||
E_coefficients, double array : stores the coefficients of the fields in the update equation for the electric field
|
||||
H_fields, double array : stores the magnetic fields of the DPW till temporal index time
|
||||
E_fields, double array : stores the electric fields of the DPW till temporal index time
|
||||
time, int : time index storing the current axis number which would be updated for the H_fields
|
||||
__________________________
|
||||
for i in range(3): # loop to update each component of the magnetic field
|
||||
materialH = G.ID[
|
||||
3 + i,
|
||||
(self.corners[0] + self.corners[3]) // 2,
|
||||
(self.corners[1] + self.corners[4]) // 2,
|
||||
(self.corners[2] + self.corners[5]) // 2,
|
||||
]
|
||||
for j in range(
|
||||
self.m[-1], self.length - self.m[-1]
|
||||
): # loop to update the magnetic field at each spatial index
|
||||
self.H_fields[i, j] = (
|
||||
G.updatecoeffsH[materialH, 0] * self.H_fields[i, j]
|
||||
+ G.updatecoeffsH[materialH, (i + 2) % 3 + 1]
|
||||
* (
|
||||
self.E_fields[(i + 1) % 3, j + self.m[(i + 2) % 3]]
|
||||
- self.E_fields[(i + 1) % 3, j]
|
||||
)
|
||||
- G.updatecoeffsH[materialH, (i + 1) % 3 + 1]
|
||||
* (
|
||||
self.E_fields[(i + 2) % 3, j + self.m[(i + 1) % 3]]
|
||||
- self.E_fields[(i + 2) % 3, j]
|
||||
)
|
||||
) # equation 8 of Tan, Potter paper
|
||||
|
||||
def update_electric_field_1D(self, G):
|
||||
"""Updates electric fields for the next time step using Equation 9 of
|
||||
DOI: 10.1109/LAWP.2009.2016851
|
||||
|
||||
Args:
|
||||
n: int stores spatial length of DPW array so that each length grid
|
||||
cell is updated when updateMagneticFields() is called.
|
||||
E_coefficients: double array stores coefficients of the fields in
|
||||
the update equation for the electric field.
|
||||
H_fields: double array stores magnetic fields of the DPW until
|
||||
temporal index time.
|
||||
E_fields: double array stores electric fields of the DPW until
|
||||
temporal index time.
|
||||
time: int time index storing current axis number which would be
|
||||
updated for the E_fields.
|
||||
|
||||
Returns:
|
||||
--------------------------
|
||||
E_fields, double array : electric field array with the axis entry for the current time added
|
||||
E_fields: double array for electric field with the axis entry for
|
||||
the current time added.
|
||||
|
||||
"""
|
||||
|
||||
for i in range(3): #loop to update each component of the electric field
|
||||
materialE = G.ID[i, (self.corners[0]+self.corners[3])//2,
|
||||
(self.corners[1]+self.corners[4])//2,
|
||||
(self.corners[2]+self.corners[5])//2 ]
|
||||
for j in range(self.m[-1], self.length-self.m[-1]): #loop to update the electric field at each spatial index
|
||||
self.E_fields[i, j] = G.updatecoeffsE[materialE, 0] * self.E_fields[i, j] + G.updatecoeffsE[materialE, (i+2)%3+1] * (
|
||||
self.H_fields[(i+2)%3, j] - self.H_fields[(i+2)%3, j-self.m[(i+1)%3]]) - G.updatecoeffsE[materialE, (i+1)%3+1] * (
|
||||
self.H_fields[(i+1)%3, j] - self.H_fields[(i+1)%3, j-self.m[(i+2)%3]]) #equation 9 of Tan, Potter paper
|
||||
|
||||
|
||||
|
||||
for i in range(3): # loop to update each component of the electric field
|
||||
materialE = G.ID[
|
||||
i,
|
||||
(self.corners[0] + self.corners[3]) // 2,
|
||||
(self.corners[1] + self.corners[4]) // 2,
|
||||
(self.corners[2] + self.corners[5]) // 2,
|
||||
]
|
||||
for j in range(
|
||||
self.m[-1], self.length - self.m[-1]
|
||||
): # loop to update the electric field at each spatial index
|
||||
self.E_fields[i, j] = (
|
||||
G.updatecoeffsE[materialE, 0] * self.E_fields[i, j]
|
||||
+ G.updatecoeffsE[materialE, (i + 2) % 3 + 1]
|
||||
* (
|
||||
self.H_fields[(i + 2) % 3, j]
|
||||
- self.H_fields[(i + 2) % 3, j - self.m[(i + 1) % 3]]
|
||||
)
|
||||
- G.updatecoeffsE[materialE, (i + 1) % 3 + 1]
|
||||
* (
|
||||
self.H_fields[(i + 1) % 3, j]
|
||||
- self.H_fields[(i + 1) % 3, j - self.m[(i + 2) % 3]]
|
||||
)
|
||||
) # equation 9 of Tan, Potter paper
|
||||
|
||||
def getField(self, i, j, k, array, m, component):
|
||||
return array[component, np.dot(m[:-1], np.array([i, j, k]))]
|
||||
|
||||
|
||||
def apply_TFSF_conditions_magnetic(self, G):
|
||||
|
||||
#**** constant x faces -- scattered-field nodes ****
|
||||
# **** constant x faces -- scattered-field nodes ****
|
||||
i = self.corners[0]
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Hy at firstX-1/2 by subtracting Ez_inc
|
||||
G.Hy[i-1, j, k] -= G.updatecoeffsH[G.ID[4, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 2)
|
||||
# correct Hy at firstX-1/2 by subtracting Ez_inc
|
||||
G.Hy[i - 1, j, k] -= G.updatecoeffsH[
|
||||
G.ID[4, i, j, k], 1
|
||||
] * self.getField(i, j, k, self.E_fields, self.m, 2)
|
||||
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Hz at firstX-1/2 by adding Ey_inc
|
||||
G.Hz[i-1, j, k] += G.updatecoeffsH[G.ID[5, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 1)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at firstX-1/2 by adding Ey_inc
|
||||
G.Hz[i - 1, j, k] += G.updatecoeffsH[
|
||||
G.ID[5, i, j, k], 1
|
||||
] * self.getField(i, j, k, self.E_fields, self.m, 1)
|
||||
|
||||
i = self.corners[3]
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Hy at lastX+1/2 by adding Ez_inc
|
||||
G.Hy[i, j, k] += G.updatecoeffsH[G.ID[4, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 2)
|
||||
# correct Hy at lastX+1/2 by adding Ez_inc
|
||||
G.Hy[i, j, k] += G.updatecoeffsH[G.ID[4, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 2
|
||||
)
|
||||
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Hz at lastX+1/2 by subtractinging Ey_inc
|
||||
G.Hz[i, j, k] -= G.updatecoeffsH[G.ID[5, i, j, k], 1] * self.getField(i, j, k, self.E_fields, self.m, 1)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at lastX+1/2 by subtractinging Ey_inc
|
||||
G.Hz[i, j, k] -= G.updatecoeffsH[G.ID[5, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 1
|
||||
)
|
||||
|
||||
#**** constant y faces -- scattered-field nodes ****
|
||||
# **** constant y faces -- scattered-field nodes ****
|
||||
j = self.corners[1]
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Hx at firstY-1/2 by adding Ez_inc
|
||||
G.Hx[i, j-1, k] += G.updatecoeffsH[G.ID[3, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 2)
|
||||
# correct Hx at firstY-1/2 by adding Ez_inc
|
||||
G.Hx[i, j - 1, k] += G.updatecoeffsH[
|
||||
G.ID[3, i, j, k], 2
|
||||
] * self.getField(i, j, k, self.E_fields, self.m, 2)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Hz at firstY-1/2 by subtracting Ex_inc
|
||||
G.Hz[i, j-1, k] -= G.updatecoeffsH[G.ID[5, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 0)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at firstY-1/2 by subtracting Ex_inc
|
||||
G.Hz[i, j - 1, k] -= G.updatecoeffsH[
|
||||
G.ID[5, i, j, k], 2
|
||||
] * self.getField(i, j, k, self.E_fields, self.m, 0)
|
||||
|
||||
j = self.corners[4]
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Hx at lastY+1/2 by subtracting Ez_inc
|
||||
G.Hx[i, j, k] -= G.updatecoeffsH[G.ID[3, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 2)
|
||||
# correct Hx at lastY+1/2 by subtracting Ez_inc
|
||||
G.Hx[i, j, k] -= G.updatecoeffsH[G.ID[3, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 2
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Hz at lastY-1/2 by adding Ex_inc
|
||||
G.Hz[i, j, k] += G.updatecoeffsH[G.ID[5, i, j, k], 2] * self.getField(i, j, k, self.E_fields, self.m, 0)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Hz at lastY-1/2 by adding Ex_inc
|
||||
G.Hz[i, j, k] += G.updatecoeffsH[G.ID[5, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 0
|
||||
)
|
||||
|
||||
#**** constant z faces -- scattered-field nodes ****
|
||||
# **** constant z faces -- scattered-field nodes ****
|
||||
k = self.corners[2]
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
#correct Hy at firstZ-1/2 by adding Ex_inc
|
||||
G.Hy[i, j, k-1] += G.updatecoeffsH[G.ID[4, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 0)
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Hy at firstZ-1/2 by adding Ex_inc
|
||||
G.Hy[i, j, k - 1] += G.updatecoeffsH[
|
||||
G.ID[4, i, j, k], 3
|
||||
] * self.getField(i, j, k, self.E_fields, self.m, 0)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
#correct Hx at firstZ-1/2 by subtracting Ey_inc
|
||||
G.Hx[i, j, k-1] -= G.updatecoeffsH[G.ID[3, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 1)
|
||||
# correct Hx at firstZ-1/2 by subtracting Ey_inc
|
||||
G.Hx[i, j, k - 1] -= G.updatecoeffsH[
|
||||
G.ID[3, i, j, k], 3
|
||||
] * self.getField(i, j, k, self.E_fields, self.m, 1)
|
||||
|
||||
k = self.corners[5]
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
#correct Hy at firstZ-1/2 by subtracting Ex_inc
|
||||
G.Hy[i, j, k] -= G.updatecoeffsH[G.ID[4, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 0)
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Hy at firstZ-1/2 by subtracting Ex_inc
|
||||
G.Hy[i, j, k] -= G.updatecoeffsH[G.ID[4, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
#correct Hx at lastZ+1/2 by adding Ey_inc
|
||||
G.Hx[i, j, k] += G.updatecoeffsH[G.ID[3, i, j, k], 3] * self.getField(i, j, k, self.E_fields, self.m, 1)
|
||||
|
||||
|
||||
|
||||
# correct Hx at lastZ+1/2 by adding Ey_inc
|
||||
G.Hx[i, j, k] += G.updatecoeffsH[G.ID[3, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.E_fields, self.m, 1
|
||||
)
|
||||
|
||||
def apply_TFSF_conditions_electric(self, G):
|
||||
#**** constant x faces -- total-field nodes ****/
|
||||
# **** constant x faces -- total-field nodes ****/
|
||||
i = self.corners[0]
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Ez at firstX face by subtracting Hy_inc
|
||||
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(i-1, j, k, self.H_fields, self.m, 1)
|
||||
# correct Ez at firstX face by subtracting Hy_inc
|
||||
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(
|
||||
i - 1, j, k, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Ey at firstX face by adding Hz_inc
|
||||
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(i-1, j, k, self.H_fields, self.m, 2)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ey at firstX face by adding Hz_inc
|
||||
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(
|
||||
i - 1, j, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
i = self.corners[3]
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Ez at lastX face by adding Hy_inc
|
||||
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(i, j, k, self.H_fields, self.m, 1)
|
||||
# correct Ez at lastX face by adding Hy_inc
|
||||
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
i = self.corners[3]
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Ey at lastX face by subtracting Hz_inc
|
||||
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(i, j, k, self.H_fields, self.m, 2)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ey at lastX face by subtracting Hz_inc
|
||||
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 1] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
#**** constant y faces -- total-field nodes ****/
|
||||
# **** constant y faces -- total-field nodes ****/
|
||||
j = self.corners[1]
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Ez at firstY face by adding Hx_inc
|
||||
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(i, j-1, k, self.H_fields, self.m, 0)
|
||||
# correct Ez at firstY face by adding Hx_inc
|
||||
G.Ez[i, j, k] += G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(
|
||||
i, j - 1, k, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Ex at firstY face by subtracting Hz_inc
|
||||
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(i, j-1, k, self.H_fields, self.m, 2)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ex at firstY face by subtracting Hz_inc
|
||||
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(
|
||||
i, j - 1, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
j = self.corners[4]
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for k in range(self.corners[2], self.corners[5]):
|
||||
#correct Ez at lastY face by subtracting Hx_inc
|
||||
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(i, j, k, self.H_fields, self.m, 0)
|
||||
# correct Ez at lastY face by subtracting Hx_inc
|
||||
G.Ez[i, j, k] -= G.updatecoeffsE[G.ID[2, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for k in range(self.corners[2], self.corners[5]+1):
|
||||
#correct Ex at lastY face by adding Hz_inc
|
||||
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(i, j, k, self.H_fields, self.m, 2)
|
||||
for k in range(self.corners[2], self.corners[5] + 1):
|
||||
# correct Ex at lastY face by adding Hz_inc
|
||||
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 2] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 2
|
||||
)
|
||||
|
||||
#**** constant z faces -- total-field nodes ****/
|
||||
# **** constant z faces -- total-field nodes ****/
|
||||
k = self.corners[2]
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
#correct Ey at firstZ face by subtracting Hx_inc
|
||||
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(i, j, k-1, self.H_fields, self.m, 0)
|
||||
# correct Ey at firstZ face by subtracting Hx_inc
|
||||
G.Ey[i, j, k] -= G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(
|
||||
i, j, k - 1, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
#correct Ex at firstZ face by adding Hy_inc
|
||||
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(i, j, k-1, self.H_fields, self.m, 1)
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Ex at firstZ face by adding Hy_inc
|
||||
G.Ex[i, j, k] += G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(
|
||||
i, j, k - 1, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
k = self.corners[5]
|
||||
for i in range(self.corners[0], self.corners[3]+1):
|
||||
for i in range(self.corners[0], self.corners[3] + 1):
|
||||
for j in range(self.corners[1], self.corners[4]):
|
||||
#correct Ey at lastZ face by adding Hx_inc
|
||||
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(i, j, k, self.H_fields, self.m, 0)
|
||||
# correct Ey at lastZ face by adding Hx_inc
|
||||
G.Ey[i, j, k] += G.updatecoeffsE[G.ID[1, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 0
|
||||
)
|
||||
|
||||
for i in range(self.corners[0], self.corners[3]):
|
||||
for j in range(self.corners[1], self.corners[4]+1):
|
||||
#correct Ex at lastZ face by subtracting Hy_inc
|
||||
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(i, j, k, self.H_fields, self.m, 1)
|
||||
|
||||
for j in range(self.corners[1], self.corners[4] + 1):
|
||||
# correct Ex at lastZ face by subtracting Hy_inc
|
||||
G.Ex[i, j, k] -= G.updatecoeffsE[G.ID[0, i, j, k], 3] * self.getField(
|
||||
i, j, k, self.H_fields, self.m, 1
|
||||
)
|
||||
|
||||
在新工单中引用
屏蔽一个用户