unable to set planar source impinging on a cylinder

[Andrea Farina]

Dear Qianquian,

I'm struggling in generating a simple configuration. A cylinder with planar illumination. I get error with mmcaddsrc, it seems that there are problems when tassellating the soruce domain. Here is my code:

%% test cylinder with planar illumination

[node, face,elem] = meshacylinder([0,0,0], [0,45,0], 10, 1);

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

elem(:,5) = 1;

%% set pattern

cfg.srctype = 'pattern';

cfg.srcparam1 = [10,0,0,4];

cfg.srcparam2 = [0,10,0,4];

cfg.srcpos = [-5, 30, 12];

cfg.srcdir = [0, 0, 1];

%% add src pattern

srcdef = struct('srctype', cfg.srctype, 'srcpos', cfg.srcpos, 'srcdir',...

    cfg.srcdir, 'srcparam1', cfg.srcparam1, 'srcparam2', cfg.srcparam2);

[node, elem] = mmcaddsrc(node, elem, ...

    mmcsrcdomain(srcdef,[min(node);max(node)], 'rotate', pi/7));

%% plot

figure(2),plotmesh(node,elem)

%% display elemprop

disp(['min(elem(:,5)) = ',num2str(min(elem(:,5)))]);

I can make it work with the last param of meshacylinder = 0.5 but then I check that no elemprop<1 are set and the second plotmesh is really strange...

Can you please help me to solve that issue? 

Thanks a lot!

Best regards

Andrea

[Fang, Qianqian]

Hi Andrea,

Sorry for the delayed reply. Had been busy and still haven't get back to you on the mmc bugs you reported on github.

For this issue, I think it was caused by the limited precision of the initial cylindrical mesh generation.

Here is the extrapolated domain that mmcaddsrc is attempting to tesselate - this is consisted of the space made by the convex-hull of the src+cylinder domain subtracting the cylinder domain.

You can debug this by setting a breakpoint at 

https://github.com/fangq/iso2mesh/blob/v1.9.8/meshrefine.m#L236

and plot with

plotmesh(allnode, bothsides);

From this plot, you see a lot of co-planar triangles overlapping each other, filling between the parallel sections along the cylinder axis direction. These "broken side panels" is because your meshacylinder() is trying to accommodate your maxvol constrain and refine the cylinder surface to break the side panel into multiple rows of triangles. Because of the finite precision, these triangles had caused qhull to incorrectly connect triangles and produces overlapping triangles.

One way to fix this, as you noticed, is to adjust the maxvol setting for meshacylinder. Use a large number, so it won't break the side panels.

I tried 100, 30, 10 or 5, all produced good results.

If you do need to refine the cylinder region to be maxvol=1, you can call meshrefine() again after the retesselation step. 

Try this

[node, face,elem] = meshacylinder([0,0,0], [0,45,0], 10, 30);

......

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

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

sizefield=(elem1(:,5) == 1)*0.3;    %% find all elements in the cylinder, set maxvol to 0.3 for each element

[node, elem] = meshrefine(node1, elem1, sizefield);

There is a possibility that such coplanar triangles can be suppressed via either convhulln()'s qhull options in matlab, or using tetgen's option, but the above method should be a workaround.

---

From: mmc-...@googlegroups.com <mmc-...@googlegroups.com> on behalf of Andrea Farina <andreaf...@gmail.com>
Sent: Thursday, October 31, 2024 4:40 PM
To: mmc-users <mmc-...@googlegroups.com>
Subject: [mmc-users] unable to set planar source impinging on a cylinder

 

--
You received this message because you are subscribed to the Google Groups "mmc-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mmc-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/mmc-users/0428470f-88bf-4c0a-ba3a-5431fd1bdf6an%40googlegroups.com.

[Andrea Farina]

Thank you very much for the help Quinquian!

Atually I noticed that changing the maxvolume to the value suggested it works only on the LINUX platform. All the values suggested fails on AppleSilicon and Windows. By putting the breakpoint at the line suggested it is clear that under Linux the plot is smooth and under Mac and Windows there are white holes and the program fails at the end of meshrefine. Here attached is an example of plot on MAC when maxvol = 30.

Is it possible that the single precision affect less the LINUX platform with respect to WIN and MAC?

Thanks a lot!

Best regards

Andrea

[Fang, Qianqian]

Hi Andrea,

I believe this behavior is also the cause of this demo_mcxyz_skinvessel bug on windows

https://github.com/fangq/mmc/issues/85

I believe this is matlab version related - if you have windows and if you use an older matlab, say 2018/2020, I think you should not see this problem. But at some point, matlab must have changed the default settings of qhull (which is the library behind convhulln used in meshrefine), and it started creating inconsistent triangulation of a lattice grid - matlab used to have this problem for all versions before R2011, somehow it is doing this again on Mac and Windows again.

Qhull has numerous options

http://www.qhull.org/html/qh-quick.htm#options

https://www.mathworks.com/help/matlab/ref/convhulln.html#mw_7d5f28dd-1dff-4374-a8e6-c4b24db2675b

I am not entirely sure what flag should be added to circumvent this problem to return to the old behavior.

---

From: mmc-...@googlegroups.com <mmc-...@googlegroups.com> on behalf of Andrea Farina <andreaf...@gmail.com>

Sent: Sunday, November 3, 2024 2:54 PM
To: mmc-users <mmc-...@googlegroups.com>
Subject: Re: [mmc-users] unable to set planar source impinging on a cylinder

 

To view this discussion visit https://groups.google.com/d/msgid/mmc-users/a509eca0-09b6-4a95-965d-62e88b7215d8n%40googlegroups.com.

[Andrea Farina]

If it can help you, before entering the convhull details, I've simply tried to compare the output of the following command on the 3 platforms:

[node, face,elem] = meshacylinder([0,0,0], [0,45,0], 10, 30);

on Linux: 

size(node) = 265  3

size(face) =    404     3

size(elem) =    927     5

on MAC and WIN

size(node) = 255  3

size(face) =    402     3

size(elem) =    864     5

Maybe is this the source of the problem?

It seems that the same command create 10 nodes less, 2 face less and 63 elements less on WIN/MAC than on LINUX...

Thanks a lot!

Regards

Andrea

[Fang, Qianqian]

Unfortunately, this is the result of the issue, not the cause.

The root cause of the problem is that the Delauney tessellation of a rectangle is ambiguous (or does not exist), because you can split the rectangle in two ways - each following one diagonal

https://www.mathworks.com/help/matlab/math/delaunay-triangulation.html

Matlab, or qhull under the hood, must set a preference to split it one way vs the other. In some version of matlab, this preference seems to be deterministic and resulting a repeatable triangulation when tessellating it multiple times. However, in the newer matlab, this heuristic seems to be lost, and it started randomly picking which way to split the square.

---

From: mmc-...@googlegroups.com <mmc-...@googlegroups.com> on behalf of Andrea Farina <andreaf...@gmail.com>

Sent: Sunday, November 3, 2024 4:06 PM

To view this discussion visit https://groups.google.com/d/msgid/mmc-users/e92107d8-8e90-4ea7-8493-b8bcb2a3b740n%40googlegroups.com.

[Andrea Farina]

Ohh this is a mess!!

Thanks a lot for reporting it!

Best regards

Andrea Farina

 

Il giorno 3 nov 2024, alle ore 22:19, 'Fang, Qianqian' via mmc-users <mmc-...@googlegroups.com> ha scritto:



To view this discussion visit https://groups.google.com/d/msgid/mmc-users/SN6PR06MB41112F9297A5E349ECEA2093E4502%40SN6PR06MB4111.namprd06.prod.outlook.com.

[Andrea Farina]

Dear Quinquian,

I've tried to follow the meshcylinder.m routine and I've found that the point where the results depend on the platform is here 

https://github.com/fangq/iso2mesh/blob/a4009b64c3220a671f3c860d8e92164925bff810/surf2mesh.m#L111

when the mex binary is called. I've checked the temporary files and this seems to be the diverging point. It seems that the tetgen.mex** file behaves differently on the different platforms.

Best regards!

Andrea

[Andrea Farina]

Dear Qinquian,

thank you again for your help and thanks a lot cite my name in the last release of MMC!

I just want to notify that with last pull from the iso2mesh repo the code you have sent me at the beginning of the threads doesn't work anymore because the meshrefine.m command gives error on tetgen. The source of error is that in the new version at lines 151 and 157 are:

moreopt = [moreopt ' -qma '];

instead of the old working 

moreopt = [moreopt ' -qa '];

It is a simple workaround, just to notify in case you want to debug.

Thanks a lot!

Kind regards

Andrea

[Qianqian Fang]

hi Andrea, sorry for the late reply - I was hoping to find time to look into this, but have been working on something else instead.

I was able to see this error, and was able to commit a fix moment ago, see my commit

https://github.com/fangq/iso2mesh/commit/68292e6b46a9d6dba41122408ea3184adeaa5229

let me know if you are able to use this patch to get around.

thanks again for reporting this issue, I would have missed!

To view this discussion visit https://groups.google.com/d/msgid/mmc-users/9a65aec0-6009-449c-9048-e558e59aeb85n%40googlegroups.com.