342 {
343 PetscFunctionBegin;
344
345 const double *elemet_size =
metaFile.ElementSize();
346 const int *dim_size =
metaFile.DimSize();
347
348 int nx = ceil(dx/ (elemet_size[0] *
sCale));
349 int ny = ceil(dy/ (elemet_size[1] *
sCale));
350 int nz = ceil(dz/ (elemet_size[2] *
sCale));
351
352
353 double x_befor_offset = x;
354 double y_befor_offset = y;
355 double z_befor_offset = z;
356
357
358
359
360 int ix = ceil(x_befor_offset/ (elemet_size[0] *
sCale));
361 int iy = ceil(y_befor_offset/ (elemet_size[1] *
sCale));
362 int iz = ceil(z_befor_offset/ (elemet_size[2] *
sCale));
363
364 vec_ix.resize(2*nx);
365 vec_iy.resize(2*ny);
366 vec_iz.resize(2*nz);
367
368 int ii = 0;
369 for(
int i = 0;
i<2*nx;
i++) {
371 if(idx >= dim_size[0] || idx < 0) continue;
372 vec_ix[ii++] = idx;
373 }
374 vec_ix.resize(ii);
375 ii = 0;
376 for(
int i = 0;
i<2*ny;
i++) {
378 if(idx >= dim_size[1] || idx < 0) continue;
379 vec_iy[ii++] = idx;
380 }
381 vec_iy.resize(ii);
382 ii = 0;
383 for(
int i = 0;
i<2*nz;
i++) {
385 if(idx >= dim_size[2] || idx < 0) continue;
386 vec_iz[ii++] = idx;
387 }
388 vec_iz.resize(ii);
389
390
391 if (dx <= 0) {
392 vec_ix.resize(1);
393 vec_ix.push_back(round(x_befor_offset/ (elemet_size[0] *
sCale)));
394 }
395 if (dy <= 0) {
396 vec_iy.resize(1);
397 vec_iy.push_back(round(y_befor_offset/ (elemet_size[1] *
sCale)));
398 }
399 if (dz <= 0) {
400 vec_iz.resize(1);
401 vec_iz.push_back(round(z_befor_offset/ (elemet_size[2] *
sCale)));
402 }
403
404
405
406 vec_idx.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
407 vec_dist.resize(vec_iz.size() * vec_ix.size() * vec_iy.size());
408 int dm_size0_dim_size1 = dim_size[0] * dim_size[1];
409
410
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430 ii = 0;
431 for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
432 double a = (*it_iz) * dm_size0_dim_size1 ;
433 double dz2 = (*it_iz*(elemet_size[2]*
sCale) - z_befor_offset) * (*it_iz*(elemet_size[2]*
sCale) - z_befor_offset);
434 for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
435 double b = (*it_iy) * dim_size[0];
436 double dy2 = (*it_iy*(elemet_size[1]*
sCale) - y_befor_offset) * (*it_iy*(elemet_size[1]*
sCale) - y_befor_offset);
437 for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
438 vec_idx[ii] =
a +
b + (*it_ix);
439 vec_dist[ii] = ( dz2 + dy2 + (*it_ix*(elemet_size[0]*
sCale) - x_befor_offset)
440 * (*it_ix*(elemet_size[0]*
sCale) - x_befor_offset) );
441 ii++;
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456 }
457 }
458 }
459
460 }
461 else {
462
463 ii = 0;
464 for(vector<int>::iterator it_iz = vec_iz.begin(); it_iz != vec_iz.end();it_iz++) {
465 double a = (*it_iz) * dm_size0_dim_size1 ;
466 for(vector<int>::iterator it_iy = vec_iy.begin(); it_iy != vec_iy.end();it_iy++) {
467 double b = (*it_iy) * dim_size[0];
468 for(vector<int>::iterator it_ix = vec_ix.begin(); it_ix != vec_ix.end();it_ix++) {
469 vec_idx[ii++] =
a +
b + (*it_ix);
470 }
471 }
472 }
473 }
474
475
476 PetscFunctionReturn(0);
477 }