search

fangq / mmc

Mesh-based Monte Carlo (MMC)
Other
40 stars 32 forks source link

Errors with tetgen when using mmcadddet #71

Closed Totrock closed 3 years ago

Totrock commented 3 years ago

Hi, when adding a second detector with mmcadddet I get error messages from tetgen.mexa64:

tetgen.mexa64: tetgen.cxx:22708: bool tetgenmesh::scoutcrossingedge(tetgenmesh::list*, tetgenmesh::list*, tetgenmesh::list*, int*): Assertion `dest(starttet) == workpt[0]' failed.
Aborted
generating tetrahedral mesh from closed surfaces ...
creating volumetric mesh from a surface mesh ...
tetgen.mexa64: tetgen.cxx:22507: void tetgenmesh::rearrangesubfaces(tetgenmesh::list*, tetgenmesh::list*, tetgenmesh::list*, int*): Assertion `matchflag == true' failed.
Aborted
error: Tetgen command failed:

I provided a minmal example in the attached files. If you change the detsize to 4 the generation works fine again. I can reproduce this on Ubuntu+Octave4.2 and Debian+Octave6.2 mmcadddet doesn't throw these errors on Windows+Matlab. Maybe it is related to: https://github.com/fangq/iso2mesh/issues/22

The mesh was generated with iso2mesh from a CT Scan of a tooth. The data is from: Prof. Karl-Heinz Kunzelmann Department of Conservative Dentistry and Periodontology, University Hospital, LMU Munich

mmcadddet_tooth_example.zip

fangq commented 3 years ago

@Totrock, thanks for reporting this. I debugged the issue, and was able to reproduce it on my Ubuntu 18.04 box.

The source of the problem is the ambiguity of creating a Delaunay triangulation of a square - basically, in the case of a square, Delaunay does not exist, all 4 nodes are on the same circle, so, there is no unique way of creating a triangulation. for example

A ------ B
|        |
|        |
C ------ D

you can create triangles "ABC and BCD" or "ABD and ACD", either way is equally fine. When you asks octave/matlab's delaunayn() function to create a Delaunay mesh from such node, you may get inconsistent results each time you call it.

this was exactly what happened - when adding your 2nd wide-field detector, the outer surface (left) and inner surface (right) of the extended space shows inconsistently tessellated surface for the square-shaped detector No.1 at the bottom.

remeshing_square_delaunay_ambiguity

The issue is not entirely octave specific. it is less frequent, but I've seen MATLAB also does this occasionally.

when your octave script failed at the last mmcadddet, you can run tetgen separately and see the errors why it failed:

cd "/tmp/iso2mesh-$USER"
fangq@/tmp/iso2mesh-fangq$ /usr/bin/tetgen1.5 -Y post_vmesh.poly
Opening post_vmesh.poly.
Delaunizing vertices...
Delaunay seconds:  0.042298
Creating surface mesh ...
Found two facets intersect each other.
  1st: [4469, 4470, 4471] #90
  2nd: [4469, 4470, 4472] #193
terminate called after throwing an instance of 'int'
Aborted (core dumped)

the triangles that bother tetgen were from that overlapping but inconsistently created triangles.

you can manually edit the .poly file to find and remove those 4 triangles, but it is not easy to do it automatically.

Totrock commented 3 years ago

Thanks for investigating and responding so quickly. I'll try your workaround.