Geometry problem definig source and detector size

[Giovanna Tramontin]

Dear Professor Fang, 

I am having a problem when defining the source and the detector with the same size as my sample, the code will not run in this case. There is some problem with the mesh retasselation.  For example, if I set my source position to [10 20 -1] and size to 24, the program works perfectly. But, it I make it bigger, but not as big as the sample (which is 50 by 50), with size 40 and position [5 5 -1], the program also does not run. The ideal situation for me would be to have the source at [0 0 -1] with size 50. 

For some cases, there is also a warning from the triangulation line: 
Warning: Some input points are not referenced by the triangulation.

This is the code: 

%
%% Geometry    
[node1, face1] = meshabox([0 0 0],[50 50 10], 2, 2);
[node2, face2] = meshacylinder([15 0+1e-1 7], [15 50-1e-1 7], 2, 2, 2);  
[node, face] = mergemesh(node1, face1, node2, face2);
[cfg.node, cfg.elem] = surf2mesh(node, face, [0 0 0], [50 50 10], 1, 2, [1,1,1;15,1,7], [], 0);

 figure(1)
 plotmesh(cfg.node, cfg.elem);

%% Source
cfg.srcpos = [0 0 -1];
cfg.srcdir = [0 0 1];  
cfg.srctype = 'pattern';
   
%Generate patterns:
val(:,:,1) = [1 2; 3 4];
M=1;

cfg.srcpattern = val;
cfg.srcparam1=[50 0 0 2];
cfg.srcparam2=[0 50 0 2];
cfg.srcnum = 1;

srcdef=struct('srctype',cfg.srctype,'srcpos',cfg.srcpos,'srcdir',cfg.srcdir,...
'srcparam1',cfg.srcparam1,'srcparam2',cfg.srcparam2);

[cfg.node,cfg.elem] = mmcaddsrc(cfg.node,cfg.elem,...
mmcsrcdomain(srcdef,[min(cfg.node);max(cfg.node)]));

figure(2)
plotmesh(cfg.node, cfg.elem);

%Detector

cfg.detpos = [0.0, 0.0, -2.0, 0.0];
cfg.detdir = [0 0 1];
cfg.detparam1 = [50.0, 0.0, 0.0, 50];
cfg.detparam2 = [0.0, 50.0, 0.0, 50];
detpos = cfg.detpos(1:3);
detdef = struct('srctype', 'planar', 'srcpos', detpos, 'srcdir', cfg.detdir, ...
         'srcparam1', cfg.detparam1, 'srcparam2', cfg.detparam2);
                   
[cfg.node, cfg.elem] = mmcadddet(cfg.node, cfg.elem, detdef);

 figure(3)
 plotmesh(cfg.node, cfg.elem);

cfg.elemprop=cfg.elem(:,5);
cfg.elem=cfg.elem(:,1:4);

cfg.e0 = tsearchn(cfg.node,cfg.elem,cfg.srcpos);

%General config
cfg.unitinmm = 1;
cfg.nphoton = 1e6;
cfg.srcpos = [10 30 -1];
cfg.srcdir = [0 0 1];
cfg.tstart = 0;                
cfg.tend = 5e-9;            
cfg.tstep = 5e-10;            
cfg.debuglevel = 'TP';
cfg.method='elem';
cfg.prop = [0, 0, 1, 1; 0.1, 0.001, 0.89, 1.37; 1, 0.01, 0.89, 1.37];

%If 1 use OpenCL, if -1, classical MMC  - THIS CHANGES THE OUTPUT
cfg.gpuid = 1;
   
%% Run simulation

cfg.issaveexit = 1; % save detected photon exit position and angles
cfg.outputtype = 'fluence';
[output, detp] = mmclab(cfg);
 

Can you help me?

Thank you for your attention, 

Giovanna

[Giovanna Tramontin]

Hello, 

I just realized that I wrote these lines two times in the code: 

cfg.srcpos = [10 30 -1];
cfg.srcdir = [0 0 1];

But still, the error, persists. 

Best wishes,

Giovanna

[Qianqian Fang]

hi Giovanna,

what operating system do you use? also on what MATLAB version?

although on my Linux machine running xubuntu 22.04 and MATLAB R2023R, I was able to produce the mesh, and launch mmclab, but it hangs at 99% - usually this suggests mesh defects or degenerated elements.

overall, I see at least two issues with this example code:

1. our mesh retessellation code was only tested on benchmarks where the wide-field source and wide-field detectors are part of the convex hull of the combined mesh. that means, they are both located on the exterior surface of the retessellated mesh. unfortunately in your use case, where the detector is on the same side as the source, and it encloses the retessellated source mesh into the interior of the convex hull (in other words, the detector is on the convex hull, but the source domain is interior). meshing such domain is tricky, because the source tetrahedron uniquely labeled (with label -1) after the mmcaddsrc() call may be further refined and split into multiple elements - which is the case in your mesh, therefore, the launch of photons may fail because the initial element no longer covers the full source aperture

2. your square-shaped pattern source aperture's circumscribing triangle has an edge that is exactly parallel to the block-shaped domain edge. this can cause problems for matlab to tessellate such domain, as a result, matlab produces nearly degenerated tetrahedra (with elem volume nearly 0)

in the below plot, you can identify such degenerated elements using

evol=elemvolume(cfg.node, cfg.elem(:,1:4));
hold on
plotmesh(cfg.node, cfg.elem(find(evol< 1e-10), :), 'facecolor','r');

I plot these problematic elements in red in the below plot overlaying on top of the final mesh after src/detector retesellations

although theses elements may be valid when processed using double-precision, in MMC, we use single-precision data to handle mesh shapes, and these degenerated elements will make ray-tracing fail (therefore, hanging mmc)

problem 2 can be potentially address by using an option 'rotate' in mmcsrcdomain(), for example

[cfg.node,cfg.elem] = mmcaddsrc(cfg.node,cfg.elem,...
    mmcsrcdomain(srcdef,[min(cfg.node);max(cfg.node)], 'rotate', pi/7));

this could avoid the parallel edges as well as degenerated elements, with the output mesh shown as below - the minimum element volume is now 4e-8.

although the above mentioned issue#2 can be mitigated using the rotation flag, both the source element (tetrahedra where photons are launched, must be labeled as -1) and detector elements (elements where escaping photons are captures. must be labeled as -2) are incorrectly labeled due to the overlapping nature of the src/det.

I think perhaps manual labeling of src and detector elements could potentially make mmclab to run, but I've never tested this. mmc does support multiple source elements (with -1 label) and dynamically compute the initial element at photon launch

https://github.com/fangq/mmc/blob/v2024.2/src/mmc_core.cl#L1305-L1342

Qianqian

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/2af54f05-962e-461c-a22e-40c5467c60fcn%40googlegroups.com.

[Giovanna Tramontin]

Dear Professor Fang, 

Thank you for your answer! 

I am using Ubuntu 20.04 and matlab '9.10.0.1684407 (R2021a) Update 3'. What happens is exactly what you described. 

I understand problem 1, but it is not the only problem, because if I put the detector and source in opposite sides, both with the same size as the sample, the code hangs at 0%. 

But for problem 2, the rotation solution solved the problem. 

Thank you for your help! 

Best wishes

Giovanna