Fwd: [mcx-users] Re: Geometry problem definig source and detector size
Qianqian Fang
forward to mmc mailing list
-------- Forwarded Message --------
| Subject: | Re: [mcx-users] Re: Geometry problem definig source and detector size |
|---|---|
| Date: | Sat, 8 Jun 2024 01:20:56 -0400 |
| From: | 'Qianqian Fang' via mcx-users <mcx-...@googlegroups.com> |
| Reply-To: | mcx-...@googlegroups.com |
| To: | mcx-...@googlegroups.com, Giovanna Tramontin <giovanna...@gmail.com> |
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
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
--Em quinta-feira, 6 de junho de 2024 às 10:19:53 UTC+2, Giovanna Tramontin escreveu:
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
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.
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/c9952a27-d057-4629-8f7f-7bec51d47152%40neu.edu.