kspaceFirstOrder3D-OMP 1.0
The C++ implementation of the k-wave toolbox for the time-domain simulation of acoustic wave fields in 3D
MatrixClasses/UXYZ_SGXYZMatrix.cpp
Go to the documentation of this file.
00001 
00033 #include <MatrixClasses/UXYZ_SGXYZMatrix.h>
00034 #include <MatrixClasses/FFTWComplexMatrix.h>
00035 
00036 
00037 using namespace std;
00038 //----------------------------------------------------------------------------//
00039 //                              Constants                                     //
00040 //----------------------------------------------------------------------------//
00041 
00042 
00043 //----------------------------------------------------------------------------//
00044 //                              Implementation                                //
00045 //                              public methods                                //
00046 //----------------------------------------------------------------------------//
00047 
00048 
00049 
00055 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2(TRealMatrix& dt_rho0_sg, TFFTWComplexMatrix& FFT){
00056     
00057     FFT.Compute_iFFT_3D_C2R(*this);
00058         
00059     
00060     const float Divider = 1.0f/(2.0f *pTotalElementCount);
00061     //dt_rho0_sgx .* real...
00062     
00063 #ifndef __NO_OMP__            
00064     #pragma omp parallel for schedule (static)     
00065 #endif
00066     for (size_t i = 0; i < pTotalElementCount; i++){
00067         pMatrixData[i] = dt_rho0_sg[i] * (pMatrixData[i] * Divider);       
00068     }
00069     
00070 
00071 }// end of Compute_dt_div_rho_sgx_mul
00072 //------------------------------------------------------------------------------
00073 
00074 
00075    
00081 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2(float dt_rho_0_sgx, TFFTWComplexMatrix& FFT){
00082     FFT.Compute_iFFT_3D_C2R(*this);
00083         
00084     
00085     const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgx;
00086     //dt_rho0_sgx .* real...
00087     
00088 #ifndef __NO_OMP__            
00089     #pragma omp parallel for schedule (static)     
00090 #endif
00091     for (size_t i = 0; i < pTotalElementCount; i++){
00092         pMatrixData[i] = pMatrixData[i] * Divider;       
00093     }
00094       
00095     
00096 }// end of Compute_dt_rho_sg_mul_ifft_div_2
00097 //------------------------------------------------------------------------------
00098    
00099 
00106 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x
00107         (float dt_rho_0_sgx, TRealMatrix & dxudxn_sgx, TFFTWComplexMatrix& FFT){
00108     
00109     
00110     
00111     FFT.Compute_iFFT_3D_C2R(*this);        
00112     
00113     const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgx;
00114     //dt_rho0_sgx .* real...
00115     
00116         
00117     #ifndef __NO_OMP__            
00118         #pragma omp for schedule (static) 
00119     #endif    
00120     for (size_t z = 0; z < pDimensionSizes.Z; z++){
00121 
00122         register size_t i = z* pDimensionSizes.Y * pDimensionSizes.X;                        
00123         for (size_t y = 0; y< pDimensionSizes.Y; y++){                                               
00124             for (size_t x = 0; x < pDimensionSizes.X; x++){              
00125                 pMatrixData[i] = pMatrixData[i] * Divider * dxudxn_sgx[x];
00126                 i++;                    
00127             } // x
00128         } // y
00129       } // z
00130 
00131     
00132     
00133     
00134 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x
00135 //------------------------------------------------------------------------------
00136 
00137 
00144 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y
00145         (float dt_rho_0_sgy, TRealMatrix & dyudyn_sgy, TFFTWComplexMatrix& FFT){
00146     
00147     
00148     
00149     FFT.Compute_iFFT_3D_C2R(*this);        
00150     
00151     const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgy;
00152     //dt_rho0_sgx .* real...
00153     
00154         
00155     #ifndef __NO_OMP__            
00156         #pragma omp for schedule (static) 
00157     #endif    
00158     for (size_t z = 0; z < pDimensionSizes.Z; z++){
00159 
00160         register size_t i = z* pDimensionSizes.Y * pDimensionSizes.X;                        
00161         for (size_t y = 0; y< pDimensionSizes.Y; y++){                                               
00162             const float dyudyn_sgy_data = dyudyn_sgy[y] * Divider;
00163             for (size_t x = 0; x < pDimensionSizes.X; x++){              
00164                 pMatrixData[i] = pMatrixData[i] * dyudyn_sgy_data;
00165                 i++;                    
00166             } // x
00167         } // y
00168       } // z
00169 
00170     
00171     
00172     
00173 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y
00174 //------------------------------------------------------------------------------
00175 
00182 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z
00183     (float dt_rho_0_sgz, TRealMatrix & dzudzn_sgz, TFFTWComplexMatrix& FFT){
00184         
00185     
00186     
00187     
00188     FFT.Compute_iFFT_3D_C2R(*this);        
00189     
00190     const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgz;
00191     //dt_rho0_sgx .* real...
00192     
00193         
00194     #ifndef __NO_OMP__            
00195         #pragma omp for schedule (static) 
00196     #endif    
00197     for (size_t z = 0; z < pDimensionSizes.Z; z++){
00198 
00199         register size_t i = z* pDimensionSizes.Y * pDimensionSizes.X;                        
00200         const float dzudzn_sgz_data = dzudzn_sgz[z] * Divider;
00201         
00202         for (size_t y = 0; y< pDimensionSizes.Y; y++){                                                           
00203             for (size_t x = 0; x < pDimensionSizes.X; x++){              
00204                 pMatrixData[i] = pMatrixData[i] * dzudzn_sgz_data;
00205                 i++;                    
00206             } // x
00207         } // y
00208       } // z
00209 
00210     
00211         
00212  }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z
00213 //------------------------------------------------------------------------------
00214    
00215 
00216 
00224 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize(TRealMatrix& FFT_p, TRealMatrix& dt_rho0, TRealMatrix& pml){
00225 
00226     const float Divider = 1.0f / pTotalElementCount;              
00227     
00228     #ifndef __NO_OMP__            
00229         #pragma omp for schedule (static) 
00230     #endif
00231     for (size_t z = 0; z < pDimensionSizes.Z; z++){        
00232         
00233         register size_t i = z* p2DDataSliceSize;        
00234         for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00235             for (size_t x = 0; x < pDimensionSizes.X; x++){
00236                                                 
00237                 register float pMatrixElement = pMatrixData[i];
00238                 
00239                 //FFT_p.ElementMultiplyMatrices(dt_rho0);
00240                 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i];
00241                 
00242                 //BSXElementRealMultiply_1D_X(abc);                
00243                 pMatrixElement *= pml[x];
00244                  
00245                  //ElementSubMatrices(FFT_p);
00246                 pMatrixElement -= FFT_p_el;
00247                 
00248                 //BSXElementRealMultiply_1D_X(abc);
00249                 pMatrixData[i] = pMatrixElement * pml[x];
00250                         
00251                 i++;
00252             } // x
00253         } // y
00254     } // z
00255     
00256 }// end of Compute_ux_sgx_normalize_Optimized
00257 //------------------------------------------------------------------------------
00258 
00265 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_uniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix& pml){
00266     
00267     const float Divider = dt_rho0 / pTotalElementCount;              
00268     
00269     #ifndef __NO_OMP__            
00270         #pragma omp for schedule (static) 
00271     #endif
00272     for (size_t z = 0; z < pDimensionSizes.Z; z++){        
00273         
00274         register size_t i = z* p2DDataSliceSize;        
00275         for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00276             for (size_t x = 0; x < pDimensionSizes.X; x++){
00277                                                 
00278                 register float pMatrixElement = pMatrixData[i];
00279                 
00280                 //FFT_p.ElementMultiplyMatrices(dt_rho0);
00281                 const float FFT_p_el = Divider * FFT_p[i];
00282                 
00283                 //BSXElementRealMultiply_1D_X(abc);                
00284                 pMatrixElement *= pml[x];
00285                  
00286                  //ElementSubMatrices(FFT_p);
00287                 pMatrixElement -= FFT_p_el;
00288                 
00289                 //BSXElementRealMultiply_1D_X(abc);
00290                 pMatrixData[i] = pMatrixElement * pml[x];
00291                         
00292                 i++;
00293             } // x
00294         } // y
00295     } // z
00296     
00297     
00298 }// end of Compute_ux_sgx_normalize
00299 //------------------------------------------------------------------------------
00300   
00301 
00309 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_nonuniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix & dxudxn_sgx, TRealMatrix& pml){
00310     
00311     const float Divider = dt_rho0 / pTotalElementCount;              
00312     
00313     #ifndef __NO_OMP__            
00314         #pragma omp for schedule (static) 
00315     #endif
00316     for (size_t z = 0; z < pDimensionSizes.Z; z++){        
00317         
00318         register size_t i = z* p2DDataSliceSize;        
00319         for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00320             for (size_t x = 0; x < pDimensionSizes.X; x++){
00321                                                 
00322                 register float pMatrixElement = pMatrixData[i];
00323                 
00324                 //FFT_p.ElementMultiplyMatrices(dt_rho0);
00325                 const float FFT_p_el = (Divider * dxudxn_sgx[x]) * FFT_p[i];
00326                 
00327                 //BSXElementRealMultiply_1D_X(abc);                
00328                 pMatrixElement *= pml[x];
00329                  
00330                  //ElementSubMatrices(FFT_p);
00331                 pMatrixElement -= FFT_p_el;
00332                 
00333                 //BSXElementRealMultiply_1D_X(abc);
00334                 pMatrixData[i] = pMatrixElement * pml[x];
00335                         
00336                 i++;
00337             } // x
00338         } // y
00339     } // z
00340     
00341     
00342 }// end of Compute_ux_sgx_normalize_scalar_nonuniform
00343 //------------------------------------------------------------------------------
00344 
00345 
00353 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize(TRealMatrix& FFT_p, TRealMatrix& dt_rho0, TRealMatrix& pml){
00354   
00355  
00356     const float Divider = 1.0f / pTotalElementCount;    
00357     #ifndef __NO_OMP__            
00358         #pragma omp for schedule (static) 
00359     #endif
00360     for (size_t z = 0; z < pDimensionSizes.Z; z++){
00361                 
00362         size_t i = z* p2DDataSliceSize;
00363         for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00364             const float pml_y = pml[y];
00365             
00366             for (size_t x = 0; x < pDimensionSizes.X; x++){
00367                 
00368                 register float pMatrixElement = pMatrixData[i];
00369                 
00370                 //FFT_p.ElementMultiplyMatrices(dt_rho0);
00371                 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i] ;
00372                 
00373                 //BSXElementRealMultiply_1D_X(abc);                
00374                 pMatrixElement *= pml_y;
00375                  
00376                  //ElementSubMatrices(FFT_p);
00377                 pMatrixElement -= FFT_p_el;
00378                 
00379                 //BSXElementRealMultiply_1D_X(abc);
00380                 pMatrixData[i] = pMatrixElement * pml_y;
00381                 
00382                 i++;
00383                         
00384             } // x
00385         } // y
00386     } // z
00387          
00388 }// end of Compute_uy_sgy_normalize_Optimized
00389 //------------------------------------------------------------------------------
00390 
00391 
00398 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_uniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix& pml){
00399     
00400     const float Divider = dt_rho0 / pTotalElementCount;    
00401     #ifndef __NO_OMP__            
00402         #pragma omp for schedule (static) 
00403     #endif
00404     for (size_t z = 0; z < pDimensionSizes.Z; z++){
00405                 
00406         size_t i = z* p2DDataSliceSize;
00407         for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00408             const float pml_y = pml[y];
00409             
00410             for (size_t x = 0; x < pDimensionSizes.X; x++){
00411                 
00412                 register float pMatrixElement = pMatrixData[i];
00413                 
00414                 //FFT_p.ElementMultiplyMatrices(dt_rho0);
00415                 const float FFT_p_el = Divider * FFT_p[i] ;
00416                 
00417                 //BSXElementRealMultiply_1D_X(abc);                
00418                 pMatrixElement *= pml_y;
00419                  
00420                  //ElementSubMatrices(FFT_p);
00421                 pMatrixElement -= FFT_p_el;
00422                 
00423                 //BSXElementRealMultiply_1D_X(abc);
00424                 pMatrixData[i] = pMatrixElement * pml_y;
00425                 
00426                 i++;
00427                         
00428             } // x
00429         } // y
00430     } // z
00431     
00432     
00433 }// end of Compute_ux_sgx_normalize
00434 //------------------------------------------------------------------------------
00435 
00436 
00444 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_nonuniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix & dyudyn_sgy, TRealMatrix& pml){
00445    
00446         const float Divider = dt_rho0 / pTotalElementCount;    
00447     #ifndef __NO_OMP__            
00448         #pragma omp for schedule (static) 
00449     #endif
00450     for (size_t z = 0; z < pDimensionSizes.Z; z++){
00451                 
00452         size_t i = z* p2DDataSliceSize;
00453         for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00454             const float pml_y = pml[y];
00455             const float dyudyn_sgy_data = dyudyn_sgy[y];
00456             for (size_t x = 0; x < pDimensionSizes.X; x++){
00457                 
00458                 register float pMatrixElement = pMatrixData[i];
00459                 
00460                 //FFT_p.ElementMultiplyMatrices(dt_rho0);
00461                 const float FFT_p_el = (Divider * dyudyn_sgy_data) * FFT_p[i] ;
00462                 
00463                 //BSXElementRealMultiply_1D_X(abc);                
00464                 pMatrixElement *= pml_y;
00465                  
00466                  //ElementSubMatrices(FFT_p);
00467                 pMatrixElement -= FFT_p_el;
00468                 
00469                 //BSXElementRealMultiply_1D_X(abc);
00470                 pMatrixData[i] = pMatrixElement * pml_y;
00471                 
00472                 i++;
00473                         
00474             } // x
00475         } // y
00476     } // z
00477     
00478     
00479 }// end of Compute_uy_sgy_normalize_scalar_nonuniform
00480 //------------------------------------------------------------------------------
00481 
00482 
00490 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize(TRealMatrix& FFT_p, TRealMatrix& dt_rho0, TRealMatrix& pml){
00491 
00492       const float Divider = 1.0f / pTotalElementCount;
00493       
00494     #ifndef __NO_OMP__            
00495         #pragma omp for schedule (static) 
00496     #endif
00497       for (size_t z = 0; z < pDimensionSizes.Z; z++){        
00498             size_t i = z* p2DDataSliceSize;
00499             const float pml_z = pml[z];
00500 
00501             for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00502                 for (size_t x = 0; x < pDimensionSizes.X; x++){                
00503 
00504                     register float pMatrixElement = pMatrixData[i];
00505 
00506                     //FFT_p.ElementMultiplyMatrices(dt_rho0);
00507                     const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i];
00508 
00509                     //BSXElementRealMultiply_1D_X(abc);                
00510                     pMatrixElement *= pml_z;
00511 
00512                      //ElementSubMatrices(FFT_p);
00513                     pMatrixElement -= FFT_p_el;
00514 
00515                     //BSXElementRealMultiply_1D_X(abc);
00516 
00517                     pMatrixData[i] = pMatrixElement * pml_z;
00518 
00519                     i++;
00520 
00521                 } // x
00522             } // y
00523         } // z
00524 
00525 }// end of Compute_uz_sgz_normalize_fast
00526 //------------------------------------------------------------------------------
00527 
00534 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_uniform(TRealMatrix& FFT_p, float& dt_rho0, TRealMatrix& pml){
00535     
00536       const float Divider = dt_rho0 / pTotalElementCount;
00537       
00538     #ifndef __NO_OMP__            
00539         #pragma omp for schedule (static) 
00540     #endif
00541       for (size_t z = 0; z < pDimensionSizes.Z; z++){        
00542             size_t i = z* p2DDataSliceSize;
00543             const float pml_z = pml[z];
00544 
00545             for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00546                 for (size_t x = 0; x < pDimensionSizes.X; x++){                
00547 
00548                     register float pMatrixElement = pMatrixData[i];
00549 
00550                     //FFT_p.ElementMultiplyMatrices(dt_rho0);
00551                     const float FFT_p_el = Divider * FFT_p[i];
00552 
00553                     //BSXElementRealMultiply_1D_X(abc);                
00554                     pMatrixElement *= pml_z;
00555 
00556                      //ElementSubMatrices(FFT_p);
00557                     pMatrixElement -= FFT_p_el;
00558 
00559                     //BSXElementRealMultiply_1D_X(abc);
00560 
00561                     pMatrixData[i] = pMatrixElement * pml_z;
00562 
00563                     i++;
00564 
00565                 } // x
00566             } // y
00567         } // z
00568 
00569     
00570     
00571 }// end of Compute_uz_sgz_normalize_scalar_uniform
00572 //------------------------------------------------------------------------------
00573 
00574 
00582 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_nonuniform(TRealMatrix& FFT_p, float& dt_rho0,TRealMatrix & dzudzn_sgz, TRealMatrix& pml){
00583      
00584     const float Divider = dt_rho0 / pTotalElementCount;
00585       
00586     #ifndef __NO_OMP__            
00587         #pragma omp for schedule (static) 
00588     #endif
00589       for (size_t z = 0; z < pDimensionSizes.Z; z++){        
00590             size_t i = z* p2DDataSliceSize;
00591             const float pml_z = pml[z];
00592             const float dzudzn_sgz_data = dzudzn_sgz[z];
00593             
00594             for (size_t y = 0; y< pDimensionSizes.Y; y++){                    
00595                 for (size_t x = 0; x < pDimensionSizes.X; x++){                
00596 
00597                     register float pMatrixElement = pMatrixData[i];
00598 
00599                     //FFT_p.ElementMultiplyMatrices(dt_rho0);
00600                     const float FFT_p_el = (Divider * dzudzn_sgz_data)* FFT_p[i];
00601 
00602                     //BSXElementRealMultiply_1D_X(abc);                
00603                     pMatrixElement *= pml_z;
00604 
00605                      //ElementSubMatrices(FFT_p);
00606                     pMatrixElement -= FFT_p_el;
00607 
00608                     //BSXElementRealMultiply_1D_X(abc);
00609 
00610                     pMatrixData[i] = pMatrixElement * pml_z;
00611 
00612                     i++;
00613 
00614                 } // x
00615             } // y
00616         } // z
00617 
00618     
00619 }// end of Compute_uz_sgz_normalize_scalar_nonuniform
00620 //------------------------------------------------------------------------------
00621 
00622 
00629 void Tuxyz_sgxyzMatrix::AddTransducerSource(TLongMatrix& u_source_index, TLongMatrix& delay_mask, TRealMatrix& transducer_signal){
00630 
00631 //#ifndef __NO_OMP__            
00632 //    #pragma omp parallel for schedule (static)     
00633 //#endif        
00634     for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++){
00635         pMatrixData[u_source_index[i]] += transducer_signal[delay_mask[i]];            
00636         delay_mask[i]++;    
00637     }
00638     
00639     
00640 }// end of AddTransducerSource
00641 //------------------------------------------------------------------------------
00642    
00643    
00644 
00645 
00656 void Tuxyz_sgxyzMatrix::Add_u_source(TRealMatrix &u_source_input, TLongMatrix & u_source_index, int t_index, long u_source_mode, long u_source_many){
00657     
00658     
00659     size_t index2D = t_index;        
00660     if (u_source_many != 0) { // is 2D
00661         index2D = t_index * u_source_index.GetTotalElementCount();
00662     }
00663     
00664     
00665     
00666     if (u_source_mode == 0){            
00667         for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++){
00668             
00669             pMatrixData[u_source_index[i]] = u_source_input[index2D];
00670             
00671             if (u_source_many != 0) index2D ++;
00672         }    
00673     }// end of dirichlet
00674     
00675     if (u_source_mode == 1){            
00676         for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++){
00677             
00678             pMatrixData[u_source_index[i]] += u_source_input[index2D];           
00679             
00680             if (u_source_many != 0) index2D ++;    
00681         }            
00682     }// end of add
00683     
00684 }// end of Add_u_source
00685 //------------------------------------------------------------------------------   
00686 
00687 
00688 
00689 
00690 
00691 
00692 
00693    
00694 //----------------------------------------------------------------------------//
00695 //                              Implementation                                //
00696 //                              private methods                                //
00697 //----------------------------------------------------------------------------//
 All Classes Files Functions Variables Typedefs Enumerations