From 94cbaabdb90662c8ab7a2061fdcbd43f16db8552 Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Tue, 19 Mar 2019 11:30:51 -0400 Subject: [PATCH] fix cylic bc, add mcxlab demo --- .../examples/demo_infinite_slab_cyclic_bc.m | 34 +++++++++++++++++++ src/mcx_core.cu | 20 ++++++----- 2 files changed, 45 insertions(+), 9 deletions(-) create mode 100644 mcxlab/examples/demo_infinite_slab_cyclic_bc.m diff --git a/mcxlab/examples/demo_infinite_slab_cyclic_bc.m b/mcxlab/examples/demo_infinite_slab_cyclic_bc.m new file mode 100644 index 00000000..fadbf63e --- /dev/null +++ b/mcxlab/examples/demo_infinite_slab_cyclic_bc.m @@ -0,0 +1,34 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang +% +% In this example, we show how to simulate an infinite slab using cyclic +% boundary condition +% +% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.space +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +clear cfg; +cfg.nphoton=1e8; +cfg.issrcfrom0=1; +cfg.vol=uint8(ones(60,60,20)); +cfg.srcdir=[0 0 1]; +cfg.gpuid=1; +cfg.autopilot=1; +cfg.prop=[0 0 1 1;0.005 1 0.8 1.37]; +cfg.tstart=0; +cfg.seed=99999; + +% a uniform planar source outside the volume +cfg.srctype='planar'; +cfg.srcpos=[0 0 0]; +cfg.srcparam1=[60 0 0 0]; +cfg.srcparam2=[0 60 0 0]; +cfg.tend=5e-9; +cfg.tstep=5e-9; +cfg.bc='ccrccr'; +flux=mcxlab(cfg); + +fcw=flux.data*cfg.tstep; +imagesc(log10(abs(squeeze(fcw(:,30,:))))) +axis equal; colorbar +title('a uniform planar source incident along an infinite slab'); diff --git a/src/mcx_core.cu b/src/mcx_core.cu index 9e24b07c..3a29c4ba 100644 --- a/src/mcx_core.cu +++ b/src/mcx_core.cu @@ -1280,15 +1280,17 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n /** launch new photon when exceed time window or moving from non-zero voxel to zero voxel without reflection */ if((mediaid==0 && (((isdet & 0xF)==0 && (!gcfg->doreflect || (gcfg->doreflect && n1==gproperty[0].w))) || (isdet==bcAbsorb || isdet==bcCylic) )) || f.t>gcfg->twin1){ if(isdet==bcCylic){ - if(flipdir==0) p.x=mcx_nextafterf(roundf(p.x+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.x: -gcfg->maxidx.x)),(v.x > 0.f)-(v.x < 0.f)); - if(flipdir==1) p.y=mcx_nextafterf(roundf(p.y+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.y: -gcfg->maxidx.y)),(v.y > 0.f)-(v.y < 0.f)); - if(flipdir==2) p.z=mcx_nextafterf(roundf(p.z+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.z: -gcfg->maxidx.z)),(v.z > 0.f)-(v.z < 0.f)); - idx1d=(int(floorf(p.z))*gcfg->dimlen.y+int(floorf(p.y))*gcfg->dimlen.x+int(floorf(p.x))); - mediaid=media[idx1d]; - isdet=mediaid & DET_MASK; /** upper 16bit is the mask of the covered detector */ - mediaid &= MED_MASK; /** lower 16bit is the medium index */ - GPUDEBUG(("Cylic boundary condition, moving photon in dir %d at %d flag, new pos=[%f %f %f]\n",flipdir,isdet,p.x,p.y,p.z)); - continue; + if(flipdir==0) p.x=mcx_nextafterf(roundf(p.x+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.x: -gcfg->maxidx.x)),(v.x > 0.f)-(v.x < 0.f)); + if(flipdir==1) p.y=mcx_nextafterf(roundf(p.y+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.y: -gcfg->maxidx.y)),(v.y > 0.f)-(v.y < 0.f)); + if(flipdir==2) p.z=mcx_nextafterf(roundf(p.z+((idx1d==OUTSIDE_VOLUME_MIN) ? gcfg->maxidx.z: -gcfg->maxidx.z)),(v.z > 0.f)-(v.z < 0.f)); + if(!(p.x<0.f||p.y<0.f||p.z<0.f||p.x>=gcfg->maxidx.x||p.y>=gcfg->maxidx.y||p.z>=gcfg->maxidx.z)){ + idx1d=(int(floorf(p.z))*gcfg->dimlen.y+int(floorf(p.y))*gcfg->dimlen.x+int(floorf(p.x))); + mediaid=media[idx1d]; + isdet=mediaid & DET_MASK; /** upper 16bit is the mask of the covered detector */ + mediaid &= MED_MASK; /** lower 16bit is the medium index */ + GPUDEBUG(("Cylic boundary condition, moving photon in dir %d at %d flag, new pos=[%f %f %f]\n",flipdir,isdet,p.x,p.y,p.z)); + continue; + } } GPUDEBUG(("direct relaunch at idx=[%d] mediaid=[%d], ref=[%d] bcflag=%d timegate=%d\n",idx1d,mediaid,gcfg->doreflect,isdet,f.t>gcfg->twin1)); if(launchnewphoton(&p,&v,&f,&rv,&prop,&idx1d,field,&mediaid,&w0,(mediaidold & DET_MASK),ppath,