fdtd_3d_pml.h

采用UMPL边界的FDTD编程

#pragma once

#include "run_enviroment.h"
#include "FDTD_1D_EH_PML_LOSS.h"

#ifndef _Struct1_
	///////////////////////////////////////////////////////////////////////////
	//Data structure definitions
	///////////////////////////////////////////////////////////////////////////
	typedef struct 
	{
		double ****A[43];   //the adresses of F G E
		long size[43][3];   //the x y and z dimesion
		char name[43][6];   //the name F(xyz)_1, G(xyz)_1, ... E(xyz)
		double **B_PML[12]; //the addresses of the PML vectors
		long size_B[12];  //the size of the PML vectors
		char name_B_PML[12][9];
		char path[128];   //the path where to save the workspace data 
						  //in case of exceeding the time limit
		int er;		      //for error handling
	}Allocate_E_comp;

	typedef struct 
	{
		double ****A[23]; //the adresses of B H
		long size[23][3]; //the x y and z dimesion
		char name[23][6]; //the name B(xyz)_1,...  H(xyz)
		double **B_PML[12]; //the addresses of the PML vectors
		long size_B[12];  //the size of the PML vectors
		char name_B_PML[12][9];
		char path[128];   //the path where to save the workspace data 
							//in case of exceeding the time limit
		int er;		      //for error handling
	}Allocate_H_comp;

	#define _Struct1_
#endif

class CFDTD_3D_PML
{
private:
	double pi, eps_0, mu_0;
	
	long ***Ind; //indicates the material type in each FDTD cell
	double **Mat; //contains the material properties  eps_r sigma
	long **Coord_ptSource; //the location of the point sources
	double **Par_ptSource; //the parameters of the point sources for real field
	
	Allocate_E_comp *Allocate_FGE_x, *Allocate_FGE_y, *Allocate_FGE_z;
	Allocate_H_comp *Allocate_BH_x, *Allocate_BH_y, *Allocate_BH_z;

	//The field components
	double ***Ex; 

	double ***Fx_1,  ***Fx_2,  ***Fx_3;
	double ***Fx_4,  ***Fx_5,  ***Fx_6;
	double ***Fx_7,  ***Fx_8,  ***Fx_9;
	double ***Fx_10, ***Fx_11, ***Fx_12;
	double ***Fx_13, ***Fx_15, ***Fx_16;
	double ***Fx_17, ***Fx_18;
	double ***Fx_19, ***Fx_20, ***Fx_21;
	double ***Fx_22, ***Fx_23, ***Fx_24;
	double ***Fx_25, ***Fx_26, ***Fx_27;
	double ***Gx_1,  ***Gx_2,  ***Gx_3,  ***Gx_7,  ***Gx_8,  ***Gx_9;
	double ***Gx_10, ***Gx_12, ***Gx_16, ***Gx_18;
	double ***Gx_19, ***Gx_20, ***Gx_21, ***Gx_25, ***Gx_26, ***Gx_27;
	
	double ***Ey;
	
	double ***Fy_1,  ***Fy_2,  ***Fy_3; 
	double ***Fy_4,  ***Fy_5,  ***Fy_6;
	double ***Fy_7,  ***Fy_8,  ***Fy_9;
	double ***Fy_10, ***Fy_11, ***Fy_12;
	double ***Fy_13, ***Fy_15, ***Fy_16;
	double ***Fy_17, ***Fy_18;
	double ***Fy_19, ***Fy_20, ***Fy_21;
	double ***Fy_22, ***Fy_23, ***Fy_24;
	double ***Fy_25, ***Fy_26, ***Fy_27;
	double ***Gy_1,  ***Gy_2,  ***Gy_3,  ***Gy_4,  ***Gy_6,  ***Gy_7,  ***Gy_8,  ***Gy_9;
	double ***Gy_19, ***Gy_20, ***Gy_21, ***Gy_22, ***Gy_24, ***Gy_25, ***Gy_26, ***Gy_27;
	
	double ***Ez;

	double ***Fz_1,  ***Fz_2,  ***Fz_3;
	double ***Fz_4,  ***Fz_5,  ***Fz_6;
	double ***Fz_7,  ***Fz_8,  ***Fz_9;
	double ***Fz_10, ***Fz_11, ***Fz_12;
	double ***Fz_13, ***Fz_15;
	double ***Fz_16, ***Fz_17, ***Fz_18;
	double ***Fz_19, ***Fz_20, ***Fz_21;
	double ***Fz_22, ***Fz_23, ***Fz_24;
	double ***Fz_25, ***Fz_26, ***Fz_27;

	double ***Gz_1,  ***Gz_3,  ***Gz_4,  ***Gz_6,  ***Gz_7,  ***Gz_9;
	double ***Gz_10, ***Gz_12, ***Gz_16, ***Gz_18;
	double ***Gz_19, ***Gz_21, ***Gz_22, ***Gz_24, ***Gz_25, ***Gz_27;

	double ***Hx;

	double ***Bx_1,  ***Bx_2,  ***Bx_3;
	double ***Bx_4,  ***Bx_6;
	double ***Bx_7,  ***Bx_8,  ***Bx_9;
	double ***Bx_10, ***Bx_12;
	double ***Bx_13, ***Bx_15;
	double ***Bx_16, ***Bx_18;
	double ***Bx_19, ***Bx_20, ***Bx_21;
	double ***Bx_22, ***Bx_24;
	double ***Bx_25, ***Bx_26, ***Bx_27;

	double ***Hy;

	double ***By_1,  ***By_2,  ***By_3;
	double ***By_4,  ***By_6;
	double ***By_7,  ***By_8,  ***By_9;
	double ***By_10, ***By_11, ***By_12;
	double ***By_16, ***By_17, ***By_18;
	double ***By_19, ***By_20, ***By_21;
	double ***By_22, ***By_24;
	double ***By_25, ***By_26, ***By_27;

	double ***Hz;

	double ***Bz_1,  ***Bz_2,  ***Bz_3;
	double ***Bz_4,  ***Bz_5,  ***Bz_6;
	double ***Bz_7,  ***Bz_8,  ***Bz_9;
	double ***Bz_10, ***Bz_12;
	double ***Bz_16, ***Bz_18;
	double ***Bz_19, ***Bz_20, ***Bz_21;
	double ***Bz_22, ***Bz_23, ***Bz_24;
	double ***Bz_25, ***Bz_26, ***Bz_27;
	
	double *K_a, *K_b;
	
	//Coefficients containing the PML boundary parameters
	//Electric field
	double *K_Gx_a_1, *K_Gx_b_1, *K_Ex_a_1, *K_Ex_b_1, *K_Ex_c_1, *K_Ex_d_1; 
	double *K_Gx_a_2, *K_Gx_b_2, *K_Ex_a_2, *K_Ex_b_2, *K_Ex_c_2, *K_Ex_d_2; 
	double *K_Gy_a_1, *K_Gy_b_1, *K_Ey_a_1, *K_Ey_b_1, *K_Ey_c_1, *K_Ey_d_1; 
	double *K_Gy_a_2, *K_Gy_b_2, *K_Ey_a_2, *K_Ey_b_2, *K_Ey_c_2, *K_Ey_d_2; 
	double *K_Gz_a_1, *K_Gz_b_1, *K_Ez_a_1, *K_Ez_b_1, *K_Ez_c_1, *K_Ez_d_1; 
	double *K_Gz_a_2, *K_Gz_b_2, *K_Ez_a_2, *K_Ez_b_2, *K_Ez_c_2, *K_Ez_d_2; 
	//Magnetic field
	double *K_Bx_a_1, *K_Bx_b_1, *K_Hx_a_1, *K_Hx_b_1, *K_Hx_c_1, *K_Hx_d_1; 
	double *K_Bx_a_2, *K_Bx_b_2, *K_Hx_a_2, *K_Hx_b_2, *K_Hx_c_2, *K_Hx_d_2; 
	double *K_By_a_1, *K_By_b_1, *K_Hy_a_1, *K_Hy_b_1, *K_Hy_c_1, *K_Hy_d_1; 
	double *K_By_a_2, *K_By_b_2, *K_Hy_a_2, *K_Hy_b_2, *K_Hy_c_2, *K_Hy_d_2; 
	double *K_Bz_a_1, *K_Bz_b_1, *K_Hz_a_1, *K_Hz_b_1, *K_Hz_c_1, *K_Hz_d_1; 
	double *K_Bz_a_2, *K_Bz_b_2, *K_Hz_a_2, *K_Hz_b_2, *K_Hz_c_2, *K_Hz_d_2; 
		
	//dimensions of the computational space
    long nx, ny, nz;
	//dimensions of the PML region
	long nPML_x_1, nPML_y_1, nPML_z_1;
	long nPML_x_2, nPML_y_2, nPML_z_2;

	//Symmetric boundaries
	int jel_Symmetric_ia, jel_Symmetric_ib, jel_Symmetric_ja;
	int jel_Symmetric_jb, jel_Symmetric_ka, jel_Symmetric_kb;
	//Anti-Symmetric boundaries
	int jel_Anti_Symmetric_ia, jel_Anti_Symmetric_ib, jel_Anti_Symmetric_ja;
	int jel_Anti_Symmetric_jb, jel_Anti_Symmetric_ka, jel_Anti_Symmetric_kb;

	//number of materials present in the investigated geometry
	long n_Mat;
	
	//number of point sources present in the computational domain
	long n_ptSource;
	int source_type;
	int jel_plane_wave;
	int pt_source_Ex, pt_source_Ey, pt_source_Ez;
	int pt_source_Hx, pt_source_Hy, pt_source_Hz;
	double constORalfa;

	//plane wave excitation with total field - scatterd field formulation
	CFDTD_1D_EH_PML_LOSS fdtd_1D;
	long n_1D, n_1D_MIN_1, nPML_1D;
	double d_1D;

	//the interface between the total field-scattered field zone
	double ct_Ex_1, ct_Ex_2, ct_Ey_1, ct_Ey_2, ct_Ez_1, ct_Ez_2;
	double ct_Hx_1, ct_Hx_2, ct_Hy_1, ct_Hy_2, ct_Hz_1, ct_Hz_2;

	//the incident field for total field scattered field formulation
	double *E_1D, *H_1D; 
	double *ll_1D_E, *ll_1D_H;

	//the direction of propagation and the polarization of the plane wave
	double teta;  //angle relative to +z axis 0 < teta < 180
	double phi;   //angle relative to +x axis 0 <= phi < 360
	double gamma; //polarization
	
	//the recuired incident magnetic fields to compute the electric field components
	double **Hz_i0, **Hy_i0, **Hz_i1, **Hy_i1;
	double **Hz_j0, **Hx_j0, **Hz_j1, **Hx_j1;
	double **Hy_k0, **Hx_k0, **Hy_k1, **Hx_k1;
	double **face_Hz_i0, **face_Hy_i0, **face_Hz_i1, **face_Hy_i1;
	double **face_Hz_j0, **face_Hx_j0, **face_Hz_j1, **face_Hx_j1;
	double **face_Hy_k0, **face_Hx_k0, **face_Hy_k1, **face_Hx_k1;

	//the recuired incident electric fields to compute the magnetic field components
	double **Ey_i0, **Ez_i0, **Ey_i1, **Ez_i1;
	double **Ex_j0, **Ez_j0, **Ex_j1, **Ez_j1;
	double **Ex_k0, **Ey_k0, **Ex_k1, **Ey_k1;
	double **face_Ey_i0, **face_Ez_i0, **face_Ey_i1, **face_Ez_i1;
	double **face_Ex_j0, **face_Ez_j0, **face_Ex_j1, **face_Ez_j1;
	double **face_Ex_k0, **face_Ey_k0, **face_Ex_k1, **face_Ey_k1;

	long n_TS; //the size of the Scattered Field zone
	int *jel_TS_planes, jel_TS;
	//the position of the Total Field - Scattered Field region 
	long n_TS_xa, n_TS_xb, n_TS_ya, n_TS_yb,n_TS_za, n_TS_zb;
	//the local position of the Total Field - Scattered Field region 
	long n_TS_xa_Loc, n_TS_xb_Loc, n_TS_ya_Loc, n_TS_yb_Loc, n_TS_za_Loc, n_TS_zb_Loc;
	long n_TS_xa_Glob, n_TS_ya_Glob, n_TS_za_Glob;
	long ista_TS,  jsta_TS,  ksta_TS; //, iend_TS, jend_TS, kend_TS;
	long l_TS_x, l_TS_y, l_TS_z;
	long l_y_Loc_HzEy_i, l_z_Loc_HzEy_i, l_y_Loc_HyEz_i, l_z_Loc_HyEz_i;
    long l_x_Loc_HzEx_j, l_z_Loc_HzEx_j, l_x_Loc_HxEz_j, l_z_Loc_HxEz_j;
	long l_x_Loc_HyEx_k, l_y_Loc_HyEx_k, l_x_Loc_HxEy_k, l_y_Loc_HxEy_k;

	//Auxiliary parameters for the Total Field - Scatterd Field formulation
	long n_TS_xa_Loc_MIN_1, n_TS_ya_Loc_MIN_1, n_TS_za_Loc_MIN_1;
	double cos_teta, sin_teta, cos_phi, sin_phi, sin_gamma, cos_gamma;
	long lyLocORlzLoc_HzEy_i, lyLocORlzLoc_HyEz_i;
	long lxLocORlzLoc_HzEx_j, lxLocORlzLoc_HxEz_j;
	long lxLocORlyLoc_HyEx_k, lxLocORlyLoc_HxEy_k;
	long l_x_Loc, l_y_Loc, l_z_Loc;
	long n_TS_xb_Loc_HzEx_j, n_TS_zb_Loc_HzEx_j, n_TS_xb_Loc_HyEx_k, n_TS_yb_Loc_HyEx_k;
	long n_TS_yb_Loc_HzEy_i, n_TS_zb_Loc_HzEy_i, n_TS_xb_Loc_HxEy_k, n_TS_yb_Loc_HxEy_k;
	long n_TS_yb_Loc_HyEz_i, n_TS_zb_Loc_HyEz_i, n_TS_xb_Loc_HxEz_j, n_TS_zb_Loc_HxEz_j;

	//time increment
	double dt;
	//elementary cell sizes
	double dx, dy, dz; 
	
	long saveFROMinst, saveTOinst;
	/////////////////////////////////////////////////////////////////////////////
	//auxiliary parameters to reduce the computational cost
	/////////////////////////////////////////////////////////////////////////////
	double TwoORpi;
	double TwoOREp0;
	double inv_TwoOREp0;
	double TwoOREp0DIVMu0;

	double inv_dx, inv_dy, inv_dz;

	double dtDIVMu0;
	double TwoOREp0DIVMu0ORdt;
	
	double inv_dxORmu_0, inv_dyORmu_0, inv_dzORmu_0;

	double dtDIVdx, dtDIVdy, dtDIVdz;
	
	double dtDIVMu0DIVdx, dtDIVMu0DIVdy, dtDIVMu0DIVdz;

	double TwoOREp0DIVMu0ORdtDIVdx, TwoOREp0DIVMu0ORdtDIVdy, TwoOREp0DIVMu0ORdtDIVdz;

	long nx_MIN_1, ny_MIN_1, nz_MIN_1;

	long nx_MIN_nPML_x_2, ny_MIN_nPML_y_2, nz_MIN_nPML_z_2;