Comment 3 for bug 898557

Apologize for the delayed answer: my university had been closed for almost the whole past week.

> how do you analyze that you have overlapping elements?

I'm attaching the generated mesh; look at tetrahedrons 8, 27948 and 22451. Do you get the same results by running the script?

In the .png, the yellow points are positioned at tetrahedron's 8 node locations.
I've stumblued upon this because dolfin /python) fails to interpolate a solution to a different (refined) mesh. And, if I understand correctly,
this is because the barycenter of tetrahedron 8 falls exactly at the boundaries of tetrahedrons 27948 and 22451. (dolfin/function/Function.cpp, Function::non_matching_eval)

A student that is working on this (he's the one who actually wrote the mesh-generation script)
claims (today) that his mesh is generated correctly; he works on a different system than mine,
so, I suspect some kind of miscompilation. Today I think I've found an issue with tetgen,
and I've written to his author. I don't know, however, if this is reason behind what I'm seeing. The e-mail is copied below.

>Maybe you can run with slightly different tetrahedron quality settings?
>
> plceditor.set_radius_edge_ratio(1.3)?
> plceditor.set_uniform_volume_constraint(somevalue)
>
>I also see that you have not included the hole on the facet with marker 3. When you do that you get another set number of nodes and elements. >Why is this out commented?

I'll look into these.

Thanks,

Marco

----------------------------
Dear Prof. Hang Si

my name is Marco Morandini.

I think I've traced down
(thanks to gcc's -Warry-bounds warning option)
a possible bug. The description below applies to tetgen 1.4.2
(the version distributed by tritetmesh,
https://launchpad.net/tritetmesh ),
but a quick check suggests that, even if line numbers are different, the same problem could still be present in version 1.4.3
(from your download site)

After the loop

    for (newtet.loc = 0; newtet.loc < 4; newtet.loc++) {

(line 22143 of tetgen.cxx)

newtet.loc will be equal to 4;

Then, with

    pointptr = org(newtet);

and

    pointptr = dest(newtet);

(lines 22185 and 22187), you will call the inline methods
triface::org and triface::dest (lines 4460 and 4464) that, in turn,
will access locver2org[t.loc][t.ver], locver2dest[t.loc][t.ver].

Unfortunately locver2org and locver2dest are

    static int locver2org[4][6];
    static int locver2dest[4][6];

Is the above analysis correct?

Thanks for your attention,

Marco Morandini