![]() |
kspaceFirstOrder3D-OMP 1.0
The C++ implementation of the k-wave toolbox for the time-domain simulation of acoustic wave fields in 3D
|
00001 00031 #include <stdexcept> 00032 00033 #include <MatrixClasses/MatrixContainer.h> 00034 00035 #include <Parameters/Parameters.h> 00036 #include <Utils/ErrorMessages.h> 00037 00038 //----------------------------------------------------------------------------// 00039 //----------------------------- CONSTANTS ------------------------------------// 00040 //----------------------------------------------------------------------------// 00041 00042 00043 00044 00045 //============================================================================// 00046 // TMatrixRecord // 00047 //============================================================================// 00048 00049 //----------------------------------------------------------------------------// 00050 //--------------------------- Public methods ---------------------------------// 00051 //----------------------------------------------------------------------------// 00052 00053 00058 TMatrixRecord::TMatrixRecord(const TMatrixRecord& src) : 00059 MatrixPtr(src.MatrixPtr), MatrixDataType(src.MatrixDataType), 00060 DimensionSizes(src.DimensionSizes), LoadData(src.LoadData), 00061 HDF5MatrixName(src.HDF5MatrixName) 00062 { 00063 00064 00065 }// end of TMatrixRecord 00066 //------------------------------------------------------------------------------ 00067 00068 00074 TMatrixRecord& TMatrixRecord::operator = (const TMatrixRecord& src){ 00075 00076 if (this != &src){ 00077 MatrixPtr = src.MatrixPtr; 00078 MatrixDataType = src.MatrixDataType; 00079 DimensionSizes = src.DimensionSizes; 00080 LoadData = src.LoadData; 00081 HDF5MatrixName = src.HDF5MatrixName; 00082 } 00083 00084 return *this; 00085 00086 }// end of operator = 00087 //------------------------------------------------------------------------------ 00088 00089 00090 00099 void TMatrixRecord::SetAllValues(TBaseMatrix * MatrixPtr, 00100 const TMatrixDataType MatrixDataType, 00101 const TDimensionSizes DimensionSizes, 00102 const bool LoadData, 00103 const string HDF5MatrixName){ 00104 00105 this->MatrixPtr = MatrixPtr; 00106 this->MatrixDataType = MatrixDataType; 00107 this->DimensionSizes = DimensionSizes; 00108 this->LoadData = LoadData; 00109 this->HDF5MatrixName = HDF5MatrixName; 00110 00111 }// end of SetAllValues 00112 //------------------------------------------------------------------------------ 00113 00114 00115 //----------------------------------------------------------------------------// 00116 //------------------------- Protected methods --------------------------------// 00117 //----------------------------------------------------------------------------// 00118 00119 //----------------------------------------------------------------------------// 00120 //-------------------------- Private methods ---------------------------------// 00121 //----------------------------------------------------------------------------// 00122 00123 00124 00125 //============================================================================// 00126 // TMatrixContainer // 00127 //============================================================================// 00128 00129 //----------------------------------------------------------------------------// 00130 //--------------------------- Public methods ---------------------------------// 00131 //----------------------------------------------------------------------------// 00132 00133 00134 00138 TMatrixContainer::~TMatrixContainer(){ 00139 00140 MatrixContainer.clear(); 00141 00142 }// end of ~TMatrixContainer 00143 //------------------------------------------------------------------------------ 00144 00145 00150 void TMatrixContainer::CreateAllObjects(){ 00151 00152 00153 for (TMatrixRecordContainer::iterator it = MatrixContainer.begin(); it != MatrixContainer.end(); it++){ 00154 00155 if (it->second.MatrixPtr != NULL) { 00156 PrintErrorAndThrowException(MatrixContainer_ERR_FMT_ReloactaionError, it->second.HDF5MatrixName, 00157 __FILE__,__LINE__); 00158 } 00159 00160 switch (it->second.MatrixDataType){ 00161 00162 case TMatrixRecord::mdtReal: { 00163 it->second.MatrixPtr = new TRealMatrix(it->second.DimensionSizes); 00164 break; 00165 } 00166 00167 case TMatrixRecord::mdtComplex: { 00168 it->second.MatrixPtr = new TComplexMatrix(it->second.DimensionSizes); 00169 break; 00170 } 00171 00172 case TMatrixRecord::mdtIndex: { 00173 it->second.MatrixPtr = new TLongMatrix(it->second.DimensionSizes); 00174 break; 00175 } 00176 case TMatrixRecord::mdtFFTW: { 00177 it->second.MatrixPtr = new TFFTWComplexMatrix(it->second.DimensionSizes); 00178 break; 00179 } 00180 00181 case TMatrixRecord::mdtUxyz: { 00182 it->second.MatrixPtr = new Tuxyz_sgxyzMatrix(it->second.DimensionSizes); 00183 break; 00184 } 00185 00186 default:{ 00187 PrintErrorAndThrowException(MatrixContainer_ERR_FMT_RecordUnknownDistributionType, it->second.HDF5MatrixName, 00188 __FILE__,__LINE__); 00189 00190 break; 00191 } 00192 00193 }// switch 00194 00195 00196 }// end for 00197 00198 }// end of CreateAllObjects 00199 //----------------------------------------------------------------------------- 00200 00201 00202 00207 void TMatrixContainer::LoadMatricesDataFromDisk(THDF5_File & HDF5_File){ 00208 00209 for (TMatrixRecordContainer::iterator it = MatrixContainer.begin(); it != MatrixContainer.end(); it++){ 00210 00211 if (it->second.LoadData) { 00212 00213 it->second.MatrixPtr->ReadDataFromHDF5File(HDF5_File, it->second.HDF5MatrixName.c_str()); 00214 00215 } 00216 } 00217 00218 }// end of LoadMatricesDataFromDisk 00219 //------------------------------------------------------------------------------ 00220 00221 00226 void TMatrixContainer::FreeAllMatrices(){ 00227 for (TMatrixRecordContainer::iterator it = MatrixContainer.begin(); it != MatrixContainer.end(); it++){ 00228 if (it->second.MatrixPtr){ 00229 delete it->second.MatrixPtr; 00230 it->second.MatrixPtr = NULL; 00231 } 00232 } 00233 00234 00235 }// end of FreeAllMatrices 00236 //------------------------------------------------------------------------------ 00237 00238 00239 00240 00245 void TMatrixContainer::AddMatricesIntoContainer(){ 00246 00247 TParameters * Params = TParameters::GetInstance(); 00248 00249 TDimensionSizes FullDim = Params->GetFullDimensionSizes(); 00250 TDimensionSizes ReducedDim = Params->GetReducedDimensionSizes(); 00251 00252 //----------------------Allocate all matrices ----------------------------// 00253 00254 MatrixContainer[kappa] .SetAllValues(NULL ,TMatrixRecord::mdtReal , ReducedDim, false, kappa_r_Name); 00255 if (!Params->Get_c0_scalar_flag()) { 00256 MatrixContainer[c2] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , true, c0_Name); 00257 } 00258 MatrixContainer[p] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, p_Name); 00259 00260 MatrixContainer[rhox] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, rhox_Name); 00261 MatrixContainer[rhoy] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, rhoy_Name); 00262 MatrixContainer[rhoz] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, rhoz_Name); 00263 00264 MatrixContainer[ux_sgx].SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, ux_sgx_Name); 00265 MatrixContainer[uy_sgy].SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, uy_sgy_Name); 00266 MatrixContainer[uz_sgz].SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, uz_sgz_Name); 00267 00268 MatrixContainer[duxdx] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, duxdx_Name); 00269 MatrixContainer[duydy] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, duydy_Name); 00270 MatrixContainer[duzdz] .SetAllValues(NULL ,TMatrixRecord::mdtReal , FullDim , false, duzdz_Name); 00271 00272 if (!Params->Get_rho0_scalar_flag()) { 00273 MatrixContainer[rho0] .SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim , true, rho0_Name); 00274 MatrixContainer[dt_rho0_sgx].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim , true, rho0_sgx_Name); 00275 MatrixContainer[dt_rho0_sgy].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim , true, rho0_sgy_Name); 00276 MatrixContainer[dt_rho0_sgz].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim , true, rho0_sgz_Name); 00277 } 00278 00279 00280 MatrixContainer[ddx_k_shift_pos].SetAllValues(NULL,TMatrixRecord::mdtComplex, TDimensionSizes(ReducedDim.X, 1, 1), true, ddx_k_shift_pos_r_Name); 00281 MatrixContainer[ddy_k_shift_pos].SetAllValues(NULL,TMatrixRecord::mdtComplex, TDimensionSizes(1, ReducedDim.Y, 1), true, ddy_k_shift_pos_Name); 00282 MatrixContainer[ddz_k_shift_pos].SetAllValues(NULL,TMatrixRecord::mdtComplex, TDimensionSizes(1, 1, ReducedDim.Z), true, ddz_k_shift_pos_Name); 00283 00284 MatrixContainer[ddx_k_shift_neg].SetAllValues(NULL,TMatrixRecord::mdtComplex, TDimensionSizes(ReducedDim.X ,1, 1), true, ddx_k_shift_neg_r_Name); 00285 MatrixContainer[ddy_k_shift_neg].SetAllValues(NULL,TMatrixRecord::mdtComplex, TDimensionSizes(1, ReducedDim.Y, 1), true, ddy_k_shift_neg_Name); 00286 MatrixContainer[ddz_k_shift_neg].SetAllValues(NULL,TMatrixRecord::mdtComplex, TDimensionSizes(1, 1, ReducedDim.Z), true, ddz_k_shift_neg_Name); 00287 00288 00289 MatrixContainer[pml_x_sgx] .SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(FullDim.X, 1, 1), true, pml_x_sgx_Name); 00290 MatrixContainer[pml_y_sgy] .SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1, FullDim.Y, 1), true, pml_y_sgy_Name); 00291 MatrixContainer[pml_z_sgz] .SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1, 1, FullDim.Z), true, pml_z_sgz_Name); 00292 00293 MatrixContainer[pml_x] .SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(FullDim.X, 1, 1), true, pml_x_Name); 00294 MatrixContainer[pml_y] .SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1, FullDim.Y, 1), true, pml_y_Name); 00295 MatrixContainer[pml_z] .SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1, 1, FullDim.Z), true, pml_z_Name); 00296 00297 if (Params->Get_nonlinear_flag()) { 00298 if (! Params->Get_BonA_scalar_flag()) { 00299 MatrixContainer[BonA] .SetAllValues (NULL,TMatrixRecord::mdtReal, FullDim , true, BonA_Name); 00300 } 00301 } 00302 00303 if (Params->Get_absorbing_flag() != 0) { 00304 if (!((Params->Get_c0_scalar_flag()) && (Params->Get_alpha_coeff_scallar_flag()))) { 00305 MatrixContainer[absorb_tau].SetAllValues (NULL,TMatrixRecord::mdtReal, FullDim , false, absorb_tau_Name); 00306 MatrixContainer[absorb_eta].SetAllValues (NULL,TMatrixRecord::mdtReal, FullDim , false, absorb_eta_Name); 00307 } 00308 MatrixContainer[absorb_nabla1].SetAllValues(NULL,TMatrixRecord::mdtReal, ReducedDim, false, absorb_nabla1_r_Name); 00309 MatrixContainer[absorb_nabla2].SetAllValues(NULL,TMatrixRecord::mdtReal, ReducedDim, false, absorb_nabla2_r_Name 00310 ); 00311 } 00312 00313 //-- index matrices --// 00314 MatrixContainer[sensor_mask_ind].SetAllValues(NULL,TMatrixRecord::mdtIndex, TDimensionSizes(1 ,1, Params->Get_sensor_mask_index_size()), true, sensor_mask_index_Name); 00315 00316 00317 00318 // if p0 source flag 00319 if (Params->Get_p0_source_flag() == 1){ 00320 MatrixContainer[p0_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim, true, p0_source_input_Name); 00321 } 00322 00323 00324 // us_index 00325 if ((Params->Get_transducer_source_flag() != 0) || 00326 (Params->Get_ux_source_flag() != 0) || 00327 (Params->Get_uy_source_flag() != 0) || 00328 (Params->Get_uz_source_flag() != 0) ){ 00329 00330 MatrixContainer[u_source_index].SetAllValues(NULL,TMatrixRecord::mdtIndex,TDimensionSizes(1 ,1, Params->Get_u_source_index_size()), true, u_source_index_Name); 00331 } 00332 00333 00334 // -- transducer source flag defined 00335 if (Params->Get_transducer_source_flag() != 0) { 00336 00337 MatrixContainer[delay_mask] .SetAllValues(NULL,TMatrixRecord::mdtIndex,TDimensionSizes(1 ,1, Params->Get_u_source_index_size()) , true, delay_mask_Name); 00338 MatrixContainer[transducer_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal ,TDimensionSizes(1 ,1, Params->Get_transducer_source_input_size()), true, transducer_source_input_Name); 00339 00340 } 00341 00342 00343 //-- p variables --// 00344 if (Params->Get_p_source_flag() != 0){ 00345 00346 if (Params->Get_p_source_many() == 0) { // 1D case 00347 00348 MatrixContainer[p_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,1, Params->Get_p_source_flag()), true, p_source_input_Name); 00349 00350 } else { // 2D case 00351 00352 MatrixContainer[p_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,Params->Get_p_source_index_size(), Params->Get_p_source_flag()), true, p_source_input_Name); 00353 } 00354 00355 00356 MatrixContainer[p_source_index].SetAllValues(NULL,TMatrixRecord::mdtIndex, TDimensionSizes(1 ,1, Params->Get_p_source_index_size()), true, p_source_index_Name); 00357 } 00358 00359 00360 00361 //----------------------------uxyz source flags---------------------------// 00362 if (Params->Get_ux_source_flag() != 0) { 00363 if (Params->Get_u_source_many() == 0) { // 1D 00364 00365 MatrixContainer[ux_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,1, Params->Get_ux_source_flag()), true, ux_source_input_Name); 00366 } 00367 else { // 2D 00368 MatrixContainer[ux_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,Params->Get_u_source_index_size(),Params->Get_ux_source_flag()), true, ux_source_input_Name); 00369 00370 } 00371 }// ux_source_input 00372 00373 00374 if (Params->Get_uy_source_flag() != 0) { 00375 if (Params->Get_u_source_many() == 0) { // 1D 00376 00377 MatrixContainer[uy_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,1, Params->Get_uy_source_flag()), true, uy_source_input_Name); 00378 } 00379 else { // 2D 00380 MatrixContainer[uy_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,Params->Get_u_source_index_size(),Params->Get_uy_source_flag()), true, uy_source_input_Name); 00381 00382 } 00383 00384 }// uy_source_input 00385 00386 if (Params->Get_uz_source_flag() != 0) { 00387 if (Params->Get_u_source_many() == 0) { // 1D 00388 00389 MatrixContainer[uz_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal,TDimensionSizes(1 ,1, Params->Get_uz_source_flag()), true, uz_source_input_Name); 00390 } 00391 else { // 2D 00392 MatrixContainer[uz_source_input].SetAllValues(NULL,TMatrixRecord::mdtReal,TDimensionSizes(1 ,Params->Get_u_source_index_size(),Params->Get_uz_source_flag()), true, uz_source_input_Name); 00393 00394 } 00395 00396 00397 }// uz_source_input 00398 00399 00400 00401 00402 //-- Non linear grid 00403 if (Params->Get_nonuniform_grid_flag()!= 0) { 00404 MatrixContainer[dxudxn].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(FullDim.X, 1, 1), true, dxudxn_Name); 00405 MatrixContainer[dyudyn].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1, FullDim.Y, 1), true, dyudyn_Name); 00406 MatrixContainer[dzudzn].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,1, FullDim.Z), true, dzudzn_Name); 00407 00408 00409 MatrixContainer[dxudxn_sgx].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(FullDim.X, 1, 1), true, dxudxn_sgx_Name); 00410 MatrixContainer[dyudyn_sgy].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1, FullDim.Y, 1), true, dyudyn_sgy_Name); 00411 MatrixContainer[dzudzn_sgz].SetAllValues(NULL,TMatrixRecord::mdtReal, TDimensionSizes(1 ,1, FullDim.Z), true, dzudzn_sgz_Name); 00412 00413 } 00414 00415 00416 //--Temporary matrices --// 00417 // this matrix used to load alpha_coeff for absorb_tau pre-calculation 00418 00419 if ((Params->Get_absorbing_flag() != 0) && (!Params->Get_alpha_coeff_scallar_flag())){ 00420 MatrixContainer[Temp_1_RS3D].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim, true, alpha_coeff_Name); 00421 }else{ 00422 MatrixContainer[Temp_1_RS3D].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim, false, ""); 00423 } 00424 00425 00426 MatrixContainer[Temp_2_RS3D].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim, false, ""); 00427 MatrixContainer[Temp_3_RS3D].SetAllValues(NULL,TMatrixRecord::mdtReal, FullDim, false, ""); 00428 00429 MatrixContainer[FFT_X_temp].SetAllValues(NULL, TMatrixRecord::mdtFFTW, ReducedDim, false, ""); 00430 MatrixContainer[FFT_Y_temp].SetAllValues(NULL, TMatrixRecord::mdtFFTW, ReducedDim, false, ""); 00431 MatrixContainer[FFT_Z_temp].SetAllValues(NULL, TMatrixRecord::mdtFFTW, ReducedDim, false, ""); 00432 00433 //--- output buffers --// 00434 00435 00436 TDimensionSizes SensorDims(Params->Get_sensor_mask_index_size(), 1, 1); 00437 00438 00439 if (Params->IsStore_p_rms()){ 00440 MatrixContainer[p_sensor_rms].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, p_rms_Name); 00441 } 00442 00443 if (Params->IsStore_p_max()){ 00444 MatrixContainer[p_sensor_max].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, p_max_Name); 00445 } 00446 00447 00448 if (Params->IsStore_u_rms()){ 00449 MatrixContainer[ux_sensor_rms].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, ux_rms_Name); 00450 MatrixContainer[uy_sensor_rms].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, uy_rms_Name); 00451 MatrixContainer[uz_sensor_rms].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, uz_rms_Name); 00452 } 00453 00454 if (Params->IsStore_u_max()){ 00455 MatrixContainer[ux_sensor_max].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, ux_max_Name); 00456 MatrixContainer[uy_sensor_max].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, uy_max_Name); 00457 MatrixContainer[uz_sensor_max].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, uz_max_Name); 00458 } 00459 00460 if (Params->IsStore_I_avg()){ 00461 MatrixContainer[Ix_sensor_avg].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, Ix_avg_Name); 00462 MatrixContainer[Iy_sensor_avg].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, Iy_avg_Name); 00463 MatrixContainer[Iz_sensor_avg].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, Iz_avg_Name); 00464 00465 } 00466 00467 if (Params->IsStore_I_max()){ 00468 MatrixContainer[Ix_sensor_max].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, Ix_max_Name); 00469 MatrixContainer[Iy_sensor_max].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, Iy_max_Name); 00470 MatrixContainer[Iz_sensor_max].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, Iz_max_Name); 00471 } 00472 00473 00474 // in case of intensity create buffer for time staggered values 00475 if ((Params->IsStore_I_avg()) || (Params->IsStore_I_max())){ 00476 MatrixContainer[p_sensor_i_1_raw] .SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, ""); 00477 MatrixContainer[ux_sensor_i_1_agr_2].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, ""); 00478 MatrixContainer[uy_sensor_i_1_agr_2].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, ""); 00479 MatrixContainer[uz_sensor_i_1_agr_2].SetAllValues(NULL,TMatrixRecord::mdtReal, SensorDims, false, ""); 00480 } 00481 00482 00483 00484 }// end of InitMatrixContainer 00485 //------------------------------------------------------------------------------ 00486 00487 00488 //----------------------------------------------------------------------------// 00489 //------------------------- Protected methods --------------------------------// 00490 //----------------------------------------------------------------------------// 00491 00492 00493 //----------------------------------------------------------------------------// 00494 //-------------------------- Private methods ---------------------------------// 00495 //----------------------------------------------------------------------------// 00496 00497 00508 void TMatrixContainer::PrintErrorAndThrowException(const char * FMT, const string HDF5MatrixName, 00509 const char * File, const int Line){ 00510 00511 fprintf(stderr,FMT, HDF5MatrixName.c_str(), File, Line); 00512 throw bad_alloc(); 00513 00514 }// end of PrintErrorAndAbort 00515 //------------------------------------------------------------------------------ 00516 00517 00518