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

- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

All Classes Files Functions Variables Typedefs Enumerations Friends Pages
UXYZ_SGXYZMatrix.cpp
Go to the documentation of this file.
1 /**
2 * @file UXYZ_SGXYZMatrix.cpp
3 * @author Jiri Jaros
4 * Faculty of Information Technology\n
5 * Brno University of Technology \n
6 * jarosjir@fit.vutbr.cz
7 *
8 * @brief The implementation file containing the particle velocity matrix.
9 *
10 * @version kspaceFirstOrder3D 2.15
11 *
12 * @date 28 July 2011, 11:37 (created) \n
13 * 26 September 2014, 14:15 (revised)
14 *
15 * @section License
16 * This file is part of the C++ extension of the k-Wave Toolbox (http://www.k-wave.org).\n
17 * Copyright (C) 2014 Jiri Jaros and Bradley Treeby
18 *
19 * This file is part of k-Wave. k-Wave is free software: you can redistribute it
20 * and/or modify it under the terms of the GNU Lesser General Public License as
21 * published by the Free Software Foundation, either version 3 of the License,
22 * or (at your option) any later version.
23 *
24 * k-Wave is distributed in the hope that it will be useful, but
25 * WITHOUT ANY WARRANTY; without even the implied warranty of
26 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
27 * See the GNU Lesser General Public License for more details.
28 *
29 * You should have received a copy of the GNU Lesser General Public License
30 * along with k-Wave. If not, see <http://www.gnu.org/licenses/>.
31 */
32
33
34 #include <MatrixClasses/UXYZ_SGXYZMatrix.h>
35 #include <MatrixClasses/FFTWComplexMatrix.h>
36
37
38 using namespace std;
39 //----------------------------------------------------------------------------//
40 // Constants //
41 //----------------------------------------------------------------------------//
42
43
44 //----------------------------------------------------------------------------//
45 // Implementation //
46 // public methods //
47 //----------------------------------------------------------------------------//
48
49
50
51 /**
52 * Compute dt ./ rho0_sgx .* the content of the matrix.
53 * @param [in] dt_rho0_sg - matrix with the component of dt .* rho0_sg{x,y,z}
54 * @param [in] FFT - FFT matrix (data will be overwritten)
55 */
56 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2(const TRealMatrix & dt_rho0_sg,
57 TFFTWComplexMatrix & FFT)
58 {
59 FFT.Compute_FFT_3D_C2R(*this);
60
61 const float Divider = 1.0f / (2.0f * pTotalElementCount);
62
63 //dt_rho0_sgx .* real...
64 #pragma omp parallel for schedule (static)
65 for (size_t i = 0; i < pTotalElementCount; i++)
66 {
67 pMatrixData[i] = dt_rho0_sg[i] * (pMatrixData[i] * Divider);
68 }
69 }// end of Compute_dt_rho_sg_mul_ifft_div_2
70 //------------------------------------------------------------------------------
71
72
73 /**
74 * Compute dt./rho0_sgx .* ifft (FFT)
75 * if rho0_sgx is scalar, uniform case
76 * @param [in] dt_rho_0_sgx - scalar value
77 * @param [in] FFT - FFT matrix (data will be overwritten)
78 */
79 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2(const float dt_rho_0_sgx,
80 TFFTWComplexMatrix& FFT)
81 {
82 FFT.Compute_FFT_3D_C2R(*this);
83
84 const float Divider = 1.0f / (2.0f * pTotalElementCount) * dt_rho_0_sgx;
85
86 //dt_rho0_sgx .* real...
87 #pragma omp parallel for schedule (static)
88 for (size_t i = 0; i < pTotalElementCount; i++)
89 {
90 pMatrixData[i] = pMatrixData[i] * Divider;
91 }
92 }// end of Compute_dt_rho_sg_mul_ifft_div_2
93 //------------------------------------------------------------------------------
94
95
96 /**
97 * Compute dt./rho0_sgx .* ifft (FFT)
98 * if rho0_sgx is scalar, nonuniform non uniform grid, x component.
99 * @param [in] dt_rho_0_sgx - scalar value
100 * @param [in] dxudxn_sgx - non-uniform mapping
101 * @param [in] FFT - FFT matrix (data will be overwritten)
102 */
103 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x
104 (const float dt_rho_0_sgx,
105 const TRealMatrix & dxudxn_sgx,
106 TFFTWComplexMatrix& FFT)
107 {
108 FFT.Compute_FFT_3D_C2R(*this);
109
110 const float Divider = 1.0f / (2.0f * pTotalElementCount) * dt_rho_0_sgx;
111
112 //dt_rho0_sgx .* real...
113 #pragma omp for schedule (static)
114 for (size_t z = 0; z < pDimensionSizes.Z; z++)
115 {
116 register size_t i = z * pDimensionSizes.Y * pDimensionSizes.X;
117 for (size_t y = 0; y < pDimensionSizes.Y; y++)
118 {
119 for (size_t x = 0; x < pDimensionSizes.X; x++)
120 {
121 pMatrixData[i] = pMatrixData[i] * Divider * dxudxn_sgx[x];
122 i++;
123 } // x
124 } // y
125 } // z
126 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x
127 //------------------------------------------------------------------------------
128
129
130 /**
131 * Compute dt./rho0_sgy .* ifft (FFT).
132 * * if rho0_sgx is scalar, nonuniform non uniform grid, y component.
133 * @param [in] dt_rho_0_sgy - scalar value
134 * @param [in] dyudyn_sgy - non-uniform mapping
135 * @param [in] FFT - FFT matrix (data will be overwritten)
136 */
137 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y
138 (const float dt_rho_0_sgy,
139 const TRealMatrix & dyudyn_sgy,
140 TFFTWComplexMatrix& FFT)
141 {
142 FFT.Compute_FFT_3D_C2R(*this);
143
144 const float Divider = 1.0f / (2.0f * pTotalElementCount) * dt_rho_0_sgy;
145
146 //dt_rho0_sgx .* real...
147 #pragma omp for schedule (static)
148 for (size_t z = 0; z < pDimensionSizes.Z; z++)
149 {
150 register size_t i = z * pDimensionSizes.Y * pDimensionSizes.X;
151 for (size_t y = 0; y < pDimensionSizes.Y; y++)
152 {
153 const float dyudyn_sgy_data = dyudyn_sgy[y] * Divider;
154 for (size_t x = 0; x < pDimensionSizes.X; x++)
155 {
156 pMatrixData[i] = pMatrixData[i] * dyudyn_sgy_data;
157 i++;
158 } // x
159 } // y
160 } // z
161 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y
162 //------------------------------------------------------------------------------
163
164 /**
165 * Compute dt./rho0_sgz .* ifft (FFT).
166 * if rho0_sgx is scalar, non uniform grid, z component.
167 * @param [in] dt_rho_0_sgz - scalar value
168 * @param [in] dzudzn_sgz - non-uniform mapping
169 * @param [in] FFT - FFT matrix (data will be overwritten)
170 */
171 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z
172 (const float dt_rho_0_sgz,
173 const TRealMatrix & dzudzn_sgz,
174 TFFTWComplexMatrix& FFT)
175 {
176 FFT.Compute_FFT_3D_C2R(*this);
177
178 const float Divider = 1.0f / (2.0f * pTotalElementCount) * dt_rho_0_sgz;
179
180 //dt_rho0_sgx .* real...
181 #pragma omp for schedule (static)
182 for (size_t z = 0; z < pDimensionSizes.Z; z++)
183 {
184 register size_t i = z * pDimensionSizes.Y * pDimensionSizes.X;
185 const float dzudzn_sgz_data = dzudzn_sgz[z] * Divider;
186
187 for (size_t y = 0; y < pDimensionSizes.Y; y++)
188 {
189 for (size_t x = 0; x < pDimensionSizes.X; x++)
190 {
191 pMatrixData[i] = pMatrixData[i] * dzudzn_sgz_data;
192 i++;
193 } // x
194 } // y
195 } // z
196 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z
197 //------------------------------------------------------------------------------
198
199
200
201 /**
202 * Compute a new value for ux_sgx.
203 * Default case (heterogenous).
204 *
205 * @param [in] FFT_p - fft of pressure
206 * @param [in] dt_rho0 - dt_rho0_sgx
207 * @param [in] pml - pml_x
208 */
209 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize(const TRealMatrix& FFT_p,
210 const TRealMatrix& dt_rho0,
211 const TRealMatrix& pml)
212 {
213 const float Divider = 1.0f / pTotalElementCount;
214
215 #pragma omp for schedule (static)
216 for (size_t z = 0; z < pDimensionSizes.Z; z++)
217 {
218 register size_t i = z * p2DDataSliceSize;
219 for (size_t y = 0; y < pDimensionSizes.Y; y++)
220 {
221 for (size_t x = 0; x < pDimensionSizes.X; x++)
222 {
223 register float pMatrixElement = pMatrixData[i];
224
225 //FFT_p.ElementMultiplyMatrices(dt_rho0);
226 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i];
227
228 //BSXElementRealMultiply_1D_X(pml);
229 pMatrixElement *= pml[x];
230
231 //ElementSubMatrices(FFT_p);
232 pMatrixElement -= FFT_p_el;
233
234 //BSXElementRealMultiply_1D_X(pml);
235 pMatrixData[i] = pMatrixElement * pml[x];
236
237 i++;
238 } // x
239 } // y
240 } // z
241 }// end of Compute_ux_sgx_normalize
242 //------------------------------------------------------------------------------
243
244 /**
245 * Compute a new value of ux_sgx.
246 * This is the case for rho0 being a scalar and a uniform grid.
247 * @param [in] FFT_p - matrix
248 * @param [in] dt_rho0 - scalar
249 * @param [in] pml - matrix
250 */
251 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_uniform(const TRealMatrix& FFT_p,
252 const float dt_rho0,
253 const TRealMatrix& pml)
254 {
255 const float Divider = dt_rho0 / pTotalElementCount;
256
257 #pragma omp for schedule (static)
258 for (size_t z = 0; z < pDimensionSizes.Z; z++)
259 {
260 register size_t i = z * p2DDataSliceSize;
261 for (size_t y = 0; y < pDimensionSizes.Y; y++)
262 {
263 for (size_t x = 0; x < pDimensionSizes.X; x++)
264 {
265 register float pMatrixElement = pMatrixData[i];
266
267 //FFT_p.ElementMultiplyMatrices(dt_rho0);
268 const float FFT_p_el = Divider * FFT_p[i];
269
270 //BSXElementRealMultiply_1D_X(pml);
271 pMatrixElement *= pml[x];
272
273 //ElementSubMatrices(FFT_p);
274 pMatrixElement -= FFT_p_el;
275
276 //BSXElementRealMultiply_1D_X(pml);
277 pMatrixData[i] = pMatrixElement * pml[x];
278
279 i++;
280 } // x
281 } // y
282 } // z
283 }// end of Compute_ux_sgx_normalize_scalar_uniform
284 //------------------------------------------------------------------------------
285
286
287 /**
288 * Compute a new value of ux_sgx.
289 * This is the case for rho0 being a scalar and a non-uniform grid.
290 * @param [in] FFT_p - matrix
291 * @param [in] dt_rho0 - scalar
292 * @param [in] dxudxn_sgx - scalar
293 * @param [in] pml - matrix
294 */
295 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_nonuniform(const TRealMatrix& FFT_p,
296 const float dt_rho0,
297 const TRealMatrix& dxudxn_sgx,
298 const TRealMatrix& pml)
299 {
300 const float Divider = dt_rho0 / pTotalElementCount;
301
302 #pragma omp for schedule (static)
303 for (size_t z = 0; z < pDimensionSizes.Z; z++)
304 {
305 register size_t i = z * p2DDataSliceSize;
306 for (size_t y = 0; y < pDimensionSizes.Y; y++)
307 {
308 for (size_t x = 0; x < pDimensionSizes.X; x++)
309 {
310 register float pMatrixElement = pMatrixData[i];
311
312 //FFT_p.ElementMultiplyMatrices(dt_rho0);
313 const float FFT_p_el = (Divider * dxudxn_sgx[x]) * FFT_p[i];
314
315 //BSXElementRealMultiply_1D_X(pml);
316 pMatrixElement *= pml[x];
317
318 //ElementSubMatrices(FFT_p);
319 pMatrixElement -= FFT_p_el;
320
321 //BSXElementRealMultiply_1D_X(pml);
322 pMatrixData[i] = pMatrixElement * pml[x];
323
324 i++;
325 } // x
326 } // y
327 } // z
328 }// end of Compute_ux_sgx_normalize_scalar_nonuniform
329 //------------------------------------------------------------------------------
330
331
332 /**
333 * Compute new value of uy_sgy.
334 * Default case (heterogenous).
335 *
336 * @param [in] FFT_p - fft of pressure
337 * @param [in] dt_rho0 - dt_rh0_sgy
338 * @param [in] pml - pml_y
339 */
340 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize(const TRealMatrix& FFT_p,
341 const TRealMatrix& dt_rho0,
342 const TRealMatrix& pml)
343 {
344 const float Divider = 1.0f / pTotalElementCount;
345
346 #pragma omp for schedule (static)
347 for (size_t z = 0; z < pDimensionSizes.Z; z++)
348 {
349 size_t i = z * p2DDataSliceSize;
350 for (size_t y = 0; y < pDimensionSizes.Y; y++)
351 {
352 const float pml_y = pml[y];
353 for (size_t x = 0; x < pDimensionSizes.X; x++)
354 {
355 register float pMatrixElement = pMatrixData[i];
356
357 //FFT_p.ElementMultiplyMatrices(dt_rho0);
358 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i];
359
360 //BSXElementRealMultiply_1D_X(pml);
361 pMatrixElement *= pml_y;
362
363 //ElementSubMatrices(FFT_p);
364 pMatrixElement -= FFT_p_el;
365
366 //BSXElementRealMultiply_1D_X(pml);
367 pMatrixData[i] = pMatrixElement * pml_y;
368
369 i++;
370 } // x
371 } // y
372 } // z
373 }// end of Compute_uy_sgy_normalize
374 //------------------------------------------------------------------------------
375
376
377 /**
378 * Compute new value of uy_sgy.
379 * This is the case for rho0 being a scalar and a uniform grid.
380 * @param [in] FFT_p - matrix
381 * @param [in] dt_rho0 - scalar
382 * @param [in] pml - matrix
383 */
384 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_uniform(const TRealMatrix& FFT_p,
385 const float dt_rho0,
386 const TRealMatrix& pml)
387 {
388 const float Divider = dt_rho0 / pTotalElementCount;
389
390 #pragma omp for schedule (static)
391 for (size_t z = 0; z < pDimensionSizes.Z; z++)
392 {
393 size_t i = z * p2DDataSliceSize;
394 for (size_t y = 0; y < pDimensionSizes.Y; y++)
395 {
396 const float pml_y = pml[y];
397 for (size_t x = 0; x < pDimensionSizes.X; x++)
398 {
399 register float pMatrixElement = pMatrixData[i];
400
401 //FFT_p.ElementMultiplyMatrices(dt_rho0);
402 const float FFT_p_el = Divider * FFT_p[i];
403
404 //BSXElementRealMultiply_1D_X(pml);
405 pMatrixElement *= pml_y;
406
407 //ElementSubMatrices(FFT_p);
408 pMatrixElement -= FFT_p_el;
409
410 //BSXElementRealMultiply_1D_X(pml);
411 pMatrixData[i] = pMatrixElement * pml_y;
412
413 i++;
414 } // x
415 } // y
416 } // z
417 }// end of Compute_uy_sgy_normalize_scalar_uniform
418 //------------------------------------------------------------------------------
419
420
421 /**
422 * Compute new value of uy_sgy,
423 * This is the case for rho0 being a scalar and a non-uniform grid.
424 * @param [in] FFT_p - matrix
425 * @param [in] dt_rho0 - scalar
426 * @param [in] dyudyn_sgy - scalar
427 * @param [in] pml - matrix
428 */
429 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_nonuniform(const TRealMatrix& FFT_p,
430 const float dt_rho0,
431 const TRealMatrix& dyudyn_sgy,
432 const TRealMatrix& pml)
433 {
434 const float Divider = dt_rho0 / pTotalElementCount;
435
436 #pragma omp for schedule (static)
437 for (size_t z = 0; z < pDimensionSizes.Z; z++)
438 {
439 size_t i = z * p2DDataSliceSize;
440 for (size_t y = 0; y < pDimensionSizes.Y; y++)
441 {
442 const float pml_y = pml[y];
443 const float dyudyn_sgy_data = dyudyn_sgy[y];
444 for (size_t x = 0; x < pDimensionSizes.X; x++)
445 {
446 register float pMatrixElement = pMatrixData[i];
447
448 //FFT_p.ElementMultiplyMatrices(dt_rho0);
449 const float FFT_p_el = (Divider * dyudyn_sgy_data) * FFT_p[i];
450
451 //BSXElementRealMultiply_1D_X(pml);
452 pMatrixElement *= pml_y;
453
454 //ElementSubMatrices(FFT_p);
455 pMatrixElement -= FFT_p_el;
456
457 //BSXElementRealMultiply_1D_X(pml);
458 pMatrixData[i] = pMatrixElement * pml_y;
459
460 i++;
461 } // x
462 } // y
463 } // z
464 }// end of Compute_uy_sgy_normalize_scalar_nonuniform
465 //------------------------------------------------------------------------------
466
467
468 /**
469 * Compute new value of uz_sgz.
470 * Default case (heterogenous).
471 *
472 * @param [in] FFT_p - fft of pressure
473 * @param [in] dt_rho0 - dt_rh0_sgz
474 * @param [in] pml - pml_z
475 */
476 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize(const TRealMatrix& FFT_p,
477 const TRealMatrix& dt_rho0,
478 const TRealMatrix& pml)
479 {
480 const float Divider = 1.0f / pTotalElementCount;
481
482 #pragma omp for schedule (static)
483 for (size_t z = 0; z < pDimensionSizes.Z; z++)
484 {
485 size_t i = z * p2DDataSliceSize;
486 const float pml_z = pml[z];
487 for (size_t y = 0; y < pDimensionSizes.Y; y++)
488 {
489 for (size_t x = 0; x < pDimensionSizes.X; x++)
490 {
491 register float pMatrixElement = pMatrixData[i];
492
493 //FFT_p.ElementMultiplyMatrices(dt_rho0);
494 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i];
495
496 //BSXElementRealMultiply_1D_X(pml);
497 pMatrixElement *= pml_z;
498
499 //ElementSubMatrices(FFT_p);
500 pMatrixElement -= FFT_p_el;
501
502 //BSXElementRealMultiply_1D_X(pml);
503 pMatrixData[i] = pMatrixElement * pml_z;
504
505 i++;
506 } // x
507 } // y
508 } // z
509 }// end of Compute_uz_sgz_normalize
510 //------------------------------------------------------------------------------
511
512 /**
513 * Compute a new value for uz_sgz.
514 * This is the case for rho0 being a scalar and a uniform grid.
515 *
516 * @param [in] FFT_p - matrix
517 * @param [in] dt_rho0 - scalar
518 * @param [in] pml - matrix
519 */
520 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_uniform(const TRealMatrix& FFT_p,
521 const float dt_rho0,
522 const TRealMatrix& pml)
523 {
524 const float Divider = dt_rho0 / pTotalElementCount;
525
526 #pragma omp for schedule (static)
527 for (size_t z = 0; z < pDimensionSizes.Z; z++)
528 {
529 size_t i = z* p2DDataSliceSize;
530 const float pml_z = pml[z];
531 for (size_t y = 0; y < pDimensionSizes.Y; y++)
532 {
533 for (size_t x = 0; x < pDimensionSizes.X; x++)
534 {
535 register float pMatrixElement = pMatrixData[i];
536
537 //FFT_p.ElementMultiplyMatrices(dt_rho0);
538 const float FFT_p_el = Divider * FFT_p[i];
539
540 //BSXElementRealMultiply_1D_X(abc);
541 pMatrixElement *= pml_z;
542
543 //ElementSubMatrices(FFT_p);
544 pMatrixElement -= FFT_p_el;
545
546 //BSXElementRealMultiply_1D_X(abc);
547
548 pMatrixData[i] = pMatrixElement * pml_z;
549
550 i++;
551 } // x
552 } // y
553 } // z
554 }// end of Compute_uz_sgz_normalize_scalar_uniform
555 //------------------------------------------------------------------------------
556
557
558 /**
559 * Compute a new value for uz_sgz.
560 * This is the case for rho0 being a scalar and a non-uniform grid.
561 * @param [in] FFT_p - matrix
562 * @param [in] dt_rho0 - scalar
563 * @param [in] dzudzn_sgz - scalar
564 * @param [in] pml - matrix
565 */
566 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_nonuniform(const TRealMatrix& FFT_p,
567 const float dt_rho0,
568 const TRealMatrix& dzudzn_sgz,
569 const TRealMatrix& pml)
570 {
571 const float Divider = dt_rho0 / pTotalElementCount;
572
573 #pragma omp for schedule (static)
574 for (size_t z = 0; z < pDimensionSizes.Z; z++)
575 {
576 size_t i = z * p2DDataSliceSize;
577 const float pml_z = pml[z];
578 const float dzudzn_sgz_data = dzudzn_sgz[z];
579
580 for (size_t y = 0; y < pDimensionSizes.Y; y++)
581 {
582 for (size_t x = 0; x < pDimensionSizes.X; x++)
583 {
584 register float pMatrixElement = pMatrixData[i];
585
586 //FFT_p.ElementMultiplyMatrices(dt_rho0);
587 const float FFT_p_el = (Divider * dzudzn_sgz_data) * FFT_p[i];
588
589 //BSXElementRealMultiply_1D_X(abc);
590 pMatrixElement *= pml_z;
591
592 //ElementSubMatrices(FFT_p);
593 pMatrixElement -= FFT_p_el;
594
595 //BSXElementRealMultiply_1D_X(abc);
596
597 pMatrixData[i] = pMatrixElement * pml_z;
598
599 i++;
600 } // x
601 } // y
602 } // z
603 }// end of Compute_uz_sgz_normalize_scalar_nonuniform
604 //------------------------------------------------------------------------------
605
606
607 /**
608 * Add transducer data to X dimension
609 * @param [in] u_source_index - Index matrix
610 * @param [in, out] delay_mask - Index matrix - all elements incremented by one
611 * @param [in] transducer_signal - Transducer signal
612 */
613 void Tuxyz_sgxyzMatrix::AddTransducerSource(const TIndexMatrix& u_source_index,
614 TIndexMatrix& delay_mask,
615 const TRealMatrix & transducer_signal)
616 {
617 #pragma omp parallel for schedule (static) if (u_source_index.GetTotalElementCount() > 1e5)
618 for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++)
619 {
620 pMatrixData[u_source_index[i]] += transducer_signal[delay_mask[i]];
621 delay_mask[i]++;
622 }
623 }// end of AddTransducerSource
624 //------------------------------------------------------------------------------
625
626
627
628
629 /**
630 * Add in velocity source terms.
631 *
632 * @param [in] u_source_input - Source input to add
633 * @param [in] u_source_index - Index matrix
634 * @param [in] t_index - Actual time step
635 * @param [in] u_source_mode - Mode 0 = dirichlet boundary, 1 = add in
636 * @param [in] u_source_many - 0 = One series, 1 = multiple series
637 *
638 */
639 void Tuxyz_sgxyzMatrix::Add_u_source(const TRealMatrix & u_source_input,
640 const TIndexMatrix& u_source_index,
641 const size_t t_index,
642 const size_t u_source_mode,
643 const size_t u_source_many)
644 {
645 size_t index2D = t_index;
646 if (u_source_many != 0)
647 { // is 2D
648 index2D = t_index * u_source_index.GetTotalElementCount();
649 }
650
651 if (u_source_mode == 0)
652 {
653 for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++)
654 {
655 pMatrixData[u_source_index[i]] = u_source_input[index2D];
656
657 if (u_source_many != 0) index2D++;
658 }
659 }// end of Dirichlet
660
661 if (u_source_mode == 1)
662 {
663 for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++)
664 {
665 pMatrixData[u_source_index[i]] += u_source_input[index2D];
666
667 if (u_source_many != 0) index2D++;
668 }
669 }// end of add
670 }// end of Add_u_source
671 //------------------------------------------------------------------------------
672
673
674 //----------------------------------------------------------------------------//
675 // Implementation //
676 // private methods //
677 //----------------------------------------------------------------------------//
Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2
void Compute_dt_rho_sg_mul_ifft_div_2(const TRealMatrix &dt_rho_0_sgx, TFFTWComplexMatrix &FFT)
Compute dt ./ rho0_sgx .* ifft (FFT).
Definition: UXYZ_SGXYZMatrix.cpp:56
Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize
void Compute_uz_sgz_normalize(const TRealMatrix &FFT_p, const TRealMatrix &dt_rho0, const TRealMatrix &pml)
Compute a new value for uz_sgz, default case.
Definition: UXYZ_SGXYZMatrix.cpp:476
Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_nonuniform
void Compute_ux_sgx_normalize_scalar_nonuniform(const TRealMatrix &FFT_p, const float dt_rho0, const TRealMatrix &dxudxn_sgx, const TRealMatrix &pml)
Compute a new value of ux_sgx, scalar, non-uniform case.
Definition: UXYZ_SGXYZMatrix.cpp:295
Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize
void Compute_uy_sgy_normalize(const TRealMatrix &FFT_p, const TRealMatrix &dt_rho0, const TRealMatrix &pml)
Compute a new value of uy_sgy, default case.
Definition: UXYZ_SGXYZMatrix.cpp:340
UXYZ_SGXYZMatrix.h
The header file containing the particle velocity matrix.
Tuxyz_sgxyzMatrix::Add_u_source
void Add_u_source(const TRealMatrix &u_source_input, const TIndexMatrix &u_source_index, const size_t t_index, const size_t u_source_mode, const size_t u_source_many)
Add in velocity source terms.
Definition: UXYZ_SGXYZMatrix.cpp:639
Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_uniform
void Compute_uz_sgz_normalize_scalar_uniform(const TRealMatrix &FFT_p, const float dt_rho0, const TRealMatrix &pml)
Compute a new value for uz_sgz, scalar, uniform case.
Definition: UXYZ_SGXYZMatrix.cpp:520
TBaseIndexMatrix::GetTotalElementCount
virtual size_t GetTotalElementCount() const
Get total element count of the matrix.
Definition: BaseIndexMatrix.h:71
Tuxyz_sgxyzMatrix::AddTransducerSource
void AddTransducerSource(const TIndexMatrix &u_source_index, TIndexMatrix &delay_mask, const TRealMatrix &transducer_signal)
Add transducer data source to X component.
Definition: UXYZ_SGXYZMatrix.cpp:613
TFFTWComplexMatrix
Class implementing 3D Real-To-Complex and Complex-To-Real transforms using FFTW interface.
Definition: FFTWComplexMatrix.h:50
Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x
void Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x(const float dt_rho_0_sgx, const TRealMatrix &dxudxn_sgx, TFFTWComplexMatrix &FFT)
Compute dt ./ rho0_sgx .* ifft (FFT), if rho0_sgx is scalar, non uniform grid, x component.
Definition: UXYZ_SGXYZMatrix.cpp:104
TFFTWComplexMatrix::Compute_FFT_3D_C2R
void Compute_FFT_3D_C2R(TRealMatrix &OutMatrix)
Compute 3D out-of-place Complex-to-Real FFT.
Definition: FFTWComplexMatrix.cpp:628
Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_nonuniform
void Compute_uy_sgy_normalize_scalar_nonuniform(const TRealMatrix &FFT_p, const float dt_rho0, const TRealMatrix &dyudyn_sgy, const TRealMatrix &pml)
Compute a new value of uy_sgy, scalar, non-uniform case.
Definition: UXYZ_SGXYZMatrix.cpp:429
Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_nonuniform
void Compute_uz_sgz_normalize_scalar_nonuniform(const TRealMatrix &FFT_p, const float dt_rho0, const TRealMatrix &dzudzn_sgz, const TRealMatrix &pml)
Compute a new value for uz_sgz, scalar, non-uniform case.
Definition: UXYZ_SGXYZMatrix.cpp:566
TRealMatrix
The class for real matrices.
Definition: RealMatrix.h:48
TIndexMatrix
The class for 64b unsigned integers (indices). It is used for sensor_mask_index or sensor_corners_mas...
Definition: IndexMatrix.h:50
Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_uniform
void Compute_uy_sgy_normalize_scalar_uniform(const TRealMatrix &FFT_p, const float dt_rho0, const TRealMatrix &pml)
Compute a new value of uy_sgy, scalar, uniform case.
Definition: UXYZ_SGXYZMatrix.cpp:384
Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z
void Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z(const float dt_rho_0_sgz, const TRealMatrix &dzudzn_sgz, TFFTWComplexMatrix &FFT)
Compute dt./rho0_sgx .* ifft (FFT), if rho0_sgx is scalar, non uniform grid, z component.
Definition: UXYZ_SGXYZMatrix.cpp:172
Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_uniform
void Compute_ux_sgx_normalize_scalar_uniform(const TRealMatrix &FFT_p, const float dt_rho0, const TRealMatrix &pml)
Compute a new value of ux_sgx, scalar, uniform case.
Definition: UXYZ_SGXYZMatrix.cpp:251
Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y
void Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y(const float dt_rho_0_sgy, const TRealMatrix &dyudyn_sgy, TFFTWComplexMatrix &FFT)
Compute dt./rho0_sgx .* ifft (FFT), if rho0_sgx is scalar, non uniform grid, y component.
Definition: UXYZ_SGXYZMatrix.cpp:138
Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize
void Compute_ux_sgx_normalize(const TRealMatrix &FFT_p, const TRealMatrix &dt_rho0, const TRealMatrix &pml)
Compute a new value of ux_sgx, default case.
Definition: UXYZ_SGXYZMatrix.cpp:209
FFTWComplexMatrix.h
The header file containing the class that implements 3D FFT using the FFTW interface.