fdtd_3d_bloch_pml.cpp

利用c++语言写的三维FDTD

		
		return 1;
	}
	
	//Reads from file the elements of the Geometry matrix 
	switch ( load_3D_Geom_long(Index, nx, ny, nz, Inp_D->path_name_Index) )
	{
		case 0:
				#ifdef USE_MPI //WIN
					cout << "process " << myrank << " : - Geom file file loaded successfull" << endl;
                #else //IBM AIX
					cout << "Geom file file loaded successfull" << endl;
				#endif	
				
				break;
		case 1:
				#ifdef USE_MPI //WIN
					cout << "process " << myrank << " : - Faild to open the Geom file" << endl;
					MPI_Finalize();
                #else //IBM AIX
					cout << "Faild to open the Geom file" << endl;
				#endif

				return 2;
		case 2:
				#ifdef USE_MPI //WIN
					cout << "process " << myrank << " : - Wrong Geometry file content" << endl;
					MPI_Finalize();
                #else //IBM AIX
					cout << "Wrong Geometry file content" << endl;
				#endif
				
				return 3;
	}
	/////////////////////////////////////////////////////////////////////////////

	double **Materials = NULL;
	long n_Mat = 2;
	Materials = Init_Matrix_2D<double>(n_Mat,3);
	if (!Materials)
	{
		#ifdef USE_MPI //WIN
			cout << "Error in the process : " << myrank << " - Memory allocation problem - **Materials" << endl;
			MPI_Finalize();
        #else//IBM AIX
			cout << "Memory allocation problem - **Materials" << endl;
		#endif
		return 1;
	}

	switch (load_2D(Materials, n_Mat, 3, Inp_D->path_name_MatParam))
	{
		case 0:
				#ifdef USE_MPI //WIN
					cout << "process " << myrank << " : Materials file loaded successfull" << endl;
                #else //IBM AIX
					cout << "Materials file loaded successfull" << endl;
				#endif	

				break;
		case 1:
				#ifdef USE_MPI //WIN
					cout << "process " << myrank << " : - Faild to open the Mat file" << endl;
					MPI_Finalize();
                #else //IBM AIX
					cout << "Faild to open the Mat file" << endl;
				#endif
				
				return 2;
		case 2:
				#ifdef USE_MPI //WIN
					cout << "process " << myrank << " : - Wrong Mat file content" << endl;
					MPI_Finalize();
                #else //IBM AIX
					cout << "Wrong Mat file content" << endl;
				#endif

				return 3;
		case 3:
				#ifdef USE_MPI //WIN
					cout << "process " << myrank << " : - Wrong Mat file dimension" << endl;
					MPI_Finalize();
                #else //IBM AIX
					cout << "Wrong Mat file dimension" << endl;
				#endif

				return 3;
	}

	double PML_eps_r_x1, PML_mu_r_x1, PML_eps_r_x2, PML_mu_r_x2;
	double PML_eps_r_y1, PML_mu_r_y1, PML_eps_r_y2, PML_mu_r_y2; 
	double PML_eps_r_z1, PML_mu_r_z1, PML_eps_r_z2, PML_mu_r_z2;

	PML_eps_r_x1 = Inp_D->PML_eps_r_x1;  PML_mu_r_x1 = Inp_D->PML_mu_r_x1;
	PML_eps_r_x2 = Inp_D->PML_eps_r_x2;  PML_mu_r_x2 = Inp_D->PML_mu_r_x2;
	PML_eps_r_y1 = Inp_D->PML_eps_r_y1;  PML_mu_r_y1 = Inp_D->PML_mu_r_y1;
	PML_eps_r_y2 = Inp_D->PML_eps_r_y2;  PML_mu_r_y2 = Inp_D->PML_mu_r_y2;
	PML_eps_r_z1 = Inp_D->PML_eps_r_z1;  PML_mu_r_z1 = Inp_D->PML_mu_r_z1;
	PML_eps_r_z2 = Inp_D->PML_eps_r_z2;  PML_mu_r_z2 = Inp_D->PML_mu_r_z2;

	/////////////////////////////////////////////////////////////////////////////
	//Initializations of the FDTD algorithm
	/////////////////////////////////////////////////////////////////////////////
	CFDTD_3D_COMPLEX FDTD_3D_compl;
	
	if(FDTD_3D_compl.Init(Index, nx, ny, nz, Materials, n_Mat, n_PML_x_1, n_PML_x_2, 
	                   n_PML_y_1, n_PML_y_2, n_PML_z_1, n_PML_z_2, dt, dx, dy, dz, 
					   PML_eps_r_x1, PML_mu_r_x1, PML_eps_r_x2, PML_mu_r_x2, 
					   PML_eps_r_y1, PML_mu_r_y1, PML_eps_r_y2, PML_mu_r_y2, 
					   PML_eps_r_z1, PML_mu_r_z1, PML_eps_r_z2, PML_mu_r_z2))
	{
		#ifdef USE_MPI //WIN
			cout << "process " << myrank << " : - Memory allocation problem - Init - FDTD_compl" << endl;
			MPI_Finalize();
        #else //IBM AIX
			cout << "Memory allocation problem - Init - FDTD_compl" << endl;
		#endif
		
		return 1;
	}

	//Init in case of symmetric or anti-symmetric boundary conditions
	FDTD_3D_compl.Init_Symmetries(Inp_D->jel_Symmetric_ia,  Inp_D->jel_Symmetric_ib,  Inp_D->jel_Symmetric_ja,
	                              Inp_D->jel_Symmetric_jb,  Inp_D->jel_Symmetric_ka,  Inp_D->jel_Symmetric_kb,
					              Inp_D->jel_Anti_Symmetric_ia, Inp_D->jel_Anti_Symmetric_ib, Inp_D->jel_Anti_Symmetric_ja,
	                              Inp_D->jel_Anti_Symmetric_jb, Inp_D->jel_Anti_Symmetric_ka, Inp_D->jel_Anti_Symmetric_kb);
	
	int *jel_TS_planes = NULL;
	if (jel_plane_wave == 0)
	{
		FDTD_3D_compl.Init_ptSource(Coord_ptSource, n_Coord, Param_Source, source_type,
			                        Inp_D->pt_source_Ex, Inp_D->pt_source_Ey, Inp_D->pt_source_Ez,
	                                Inp_D->pt_source_Hx, Inp_D->pt_source_Hy, Inp_D->pt_source_Hz);
	}
	else
	{
		jel_TS_planes = (int *) calloc(6,sizeof(int));
		if (jel_TS_planes)
		{
			#ifdef USE_MPI //WIN
				cout << "process " << myrank << " : - Memory allocation problem - Init - FDTD_compl" << endl;
				MPI_Finalize();
			#else //IBM AIX
				cout << "Memory allocation problem - Init - FDTD_compl" << endl;
			#endif
		}
		
		if (Inp_D->jel_TS_xa == 1)
		{
			jel_TS_planes[0] = 1;
		}

		if (Inp_D->jel_TS_xb == 1)
		{
			jel_TS_planes[1] = 1;
		}

		if (Inp_D->jel_TS_ya == 1)
		{
			jel_TS_planes[2] = 1;
		}

		if (Inp_D->jel_TS_yb == 1)
		{
			jel_TS_planes[3] = 1;
		}

		if (Inp_D->jel_TS_za == 1)
		{
			jel_TS_planes[4] = 1;
		}

		if (Inp_D->jel_TS_zb == 1)
		{
			jel_TS_planes[5] = 1;
		}

		FDTD_3D_compl.Init_TS(teta, phi, gamma, n_TS, Inp_D->nPML_1D, source_type, X0, t0,
		                      tw, omega, phase, Inp_D->PML_eps_r_1D, Inp_D->PML_mu_r_1D, 
							  jel_TS_planes);

	}
	
	/////////////////////////////////////////////////////////////////////////////
	//Initializations for data save
	/////////////////////////////////////////////////////////////////////////////

	//Create the directory to save output data
	char *Path = NULL, *Path_Data = NULL;
	Path =(char *) calloc(512,sizeof(char));
	Path_Data =(char *) calloc(512,sizeof(char));
	if (!Path || !Path_Data)
	{
		#ifdef USE_MPI //WIN
			cout << "process " << myrank << " : - Memory allocation problem - *Path or *Path_Data" << endl;
			MPI_Finalize();
        #else //IBM AIX
			cout << "Memory allocation problem - *Path or *Path_Data" << endl;
		#endif

		return 1;
	}

	#ifdef run_enviroment 
		//for WIN
		GetCurrentDirectory(512, Path); 
		strcpy(Path_Data,Path);
		strcat(Path_Data,"/data");
		CreateDirectory(Path_Data,0);
	#else
		//for UNIX-AIX
		int n_Path_Data = 0;
		long cont;
		mode_t Mode = S_IRWXU;
		if (myrank == 0)
		{
			strcpy(Path,Inp_D->path_save_data);
			strcpy(Path_Data,Path);
			cont = 1;
			while (mkdir(Path_Data, Mode)) //create the main directory
			{
				if(errno == EEXIST) //check if the directory name exists 
				{
					sprintf(Path_Data, "%s_%i", Path, cont);
				}
				else
				{
					cout << "Directory creation problem: " << Path_Data << endl;
					MPI_Finalize();
					return 1;
				}
				cont++;
			}
			
			n_Path_Data = strlen(Path_Data);
			
			cout << Path_Data << endl;
		}
		MPI_Bcast(&n_Path_Data,1,MPI_INT,0,MPI_COMM_WORLD);
		MPI_Bcast(Path_Data,n_Path_Data,MPI_CHAR,0,MPI_COMM_WORLD);
	#endif
	
	int save_field_slices = Inp_D->save_field; //to save data - Ex, Ey, Ez, Hx, Hy, Hz slices
	if (save_field_slices == 1)
	{
		#ifdef run_enviroment 
			//for WIN
			sprintf(Path, "%s/data_EH_re", Path_Data);
			CreateDirectory(Path,0);//for WIN

			sprintf(Path, "%s/data_EH_re/data_Ex", Path_Data);
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_re/data_Ey");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_re/data_Ez");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_re/data_Hx");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_re/data_Hy");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_re/data_Hz");
			CreateDirectory(Path,0);//for WIN

			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_im");
			CreateDirectory(Path,0);//for WIN

			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_im/data_Ex");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_im/data_Ey");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_im/data_Ez");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_im/data_Hx");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_im/data_Hy");
			CreateDirectory(Path,0);//for WIN
			
			strcpy(Path,Path_Data);
			strcat(Path,"/data_EH_im/data_Hz");
			CreateDirectory(Path,0);//for WIN

			sprintf(Path, "%s", Path_Data);
		#else
			//for UNIX-AIX
			sprintf(Path, "%s/%i", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX

			sprintf(Path, "%s/%i/data_EH_re", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX

			sprintf(Path, "%s/%i/data_EH_re/data_Ex", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
									
			sprintf(Path, "%s/%i/data_EH_re/data_Ey", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_re/data_Ez", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_re/data_Hx", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_re/data_Hy", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_re/data_Hz", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX

			sprintf(Path, "%s/%i/data_EH_im", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX

			sprintf(Path, "%s/%i/data_EH_im/data_Ex", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
									
			sprintf(Path, "%s/%i/data_EH_im/data_Ey", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_im/data_Ez", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_im/data_Hx", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_im/data_Hy", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
						
			sprintf(Path, "%s/%i/data_EH_im/data_Hz", Path_Data, myrank);
			mkdir(Path,Mode);//for UNIX-AIX
	
			sprintf(Path, "%s/%i", Path_Data, myrank);
		#endif
		//middle planes of the computational space - for data save
		long nx_yz = (long) (nx/2);
		long ny_xz = (long) (ny/2); 
		long nz_xy = (long) (nz/2);	

		long n_x_a = Inp_D->n_x_a;
		long n_x_b = Inp_D->n_x_b;
		long n_y_a = Inp_D->n_y_a;
		long n_y_b = Inp_D->n_y_b;
		long n_z_a = Inp_D->n_z_a;
		long n_z_b = Inp_D->n_z_b;

		long nr_Save = 10;

		long saveFROMinst = Inp_D->saveFROMinst;
		long saveTOinst   = Inp_D->saveTOinst;

		if (FDTD_3D_compl.Init_Save_FieldSlices(n_x_a, n_x_b, n_y_a, n_y_b, n_z_a, n_z_b, 
				 						        nz_xy, ny_xz, nx_yz, nr_Save, saveFROMinst,
												saveTOinst, Path))
		{
			#ifdef USE_MPI //WIN
				cout << "process " << myrank << " : - Memory allocation problem - Init_Save_FieldSlices - FDTD" << endl;
				MPI_Finalize();
			#else //IBM AIX
				cout << "Memory allocation problem - Init_Save_FieldSlices - FDTD" << endl;
			#endif
			
			return 1;
		}
	} 
	//Set the points for the followed field components
	long n_Ind_F;
	long **Ind_F = NULL;
	if (Inp_D->load_foll_field_points == 1)
	{
		n_Ind_F = Inp_D->n_Ind_F;
		Ind_F = Init_Matrix_2D<long >(n_Ind_F,3);
		if(!Ind_F)
		{
			#ifdef USE_MPI //WIN
				cout << "process " << myrank << " : - Memory allocation problem - Ind_F" << endl;
				MPI_Finalize();
			#else //IBM AIX
				cout << "Memory allocation problem - Ind_F" << endl;
			#endif