| kspaceFirstOrder3D-OMP 1.0 The C++ implementation of the k-wave toolbox for the time-domain simulation of acoustic wave fields in 3D

- Main Page
- Classes
- Files

- File List
- File Members

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