libmesh icon indicating copy to clipboard operation
libmesh copied to clipboard

fix the quadratic element issue in TriangulatorInterface

Open miaoyinb opened this issue 10 months ago • 11 comments

closes #4144

miaoyinb avatar Apr 18 '25 17:04 miaoyinb

The original approach uses node ID as the key to find the corresponding midpoint. The ID seems to have changed during triangulation.

miaoyinb avatar Apr 18 '25 17:04 miaoyinb

Job Coverage, step Generate coverage on 33866f9 wanted to post the following:

Coverage

e50c9b #4145 33866f
Total Total +/- New
Rate 63.42% 63.42% +0.00% 100.00%
Hits 74838 74842 +4 2
Misses 43166 43163 -3 0

Diff coverage report

Full coverage report

This comment will be updated on new commits.

moosebuild avatar Apr 18 '25 20:04 moosebuild

This invalidates the comment explaining all_midpoints:

  // [...] all_midpoints[{p, m}] will be the mth midpoint
  // location following after point p (when traversing a triangle
  // counter-clockwise)

And it doesn't change uses of all_midpoints that clearly make that assumption too, explicitly querying elem->point(n) for n<3 to get the point for the key. None of them needed to be changed?

I think we need a new unit test to demonstrate the bug and cover the fix.

The only change this PR makes is to store the node point coordinates instead of the node id of the first node when a segment (i.e. two nodes) is created. Later on, this point coordinates are used as p in all_midpoints[{p, m}] instead of using _mesh.point() and the stored node id. Thus, I do not think anything else needs to change; the only issue of the original code is that the node ids change during triangulation.

I agree that a unit test would help.

miaoyinb avatar Apr 22 '25 03:04 miaoyinb

Ah, I get it now. I saw segment_midpoints_keys and jumped to the conclusion that we were storing midpoints, rather than non-midpoint-points-used-as-keys-to-find-midpoints. We've got to get a clearer name for that, at the very least, and clearer comments explaining what we're doing and why.

To do that, of course, we have to understand why. I don't like fixing the effects of a problem whose cause I don't understand. It's like putting a bucket to catch water under a drip from the ceiling. Your carpet is no longer getting wet, that's great, but the drip shouldn't have been there and it might be damaging your drywall if you consider the problem solved instead of investigating further.

We're creating this->segments in elems_to_segments(), at the same time we're giving them segment_midpoints. How are we losing the connection between them afterward? Is this an insert_any_extra_boundary_points() case? Are we calling a prepare_for_use() somewhere I don't see that's renumbering nodes in the middle? We use segments repeatedly in insert_refinement_points() too, so if something with it is getting screwed up on iteration N then I don't want to just go on to N+1 because we've worked around the use case in increase_triangle_order().

roystgnr avatar Apr 23 '25 14:04 roystgnr

Could you come up with a unit test? Aside from making sure we don't regress on this, that should be enough for me to go through in the debugger and find where and why we have node ids changing during triangulation.

roystgnr avatar Apr 23 '25 14:04 roystgnr

Sure, I'll aim to work on a unit test later this week, once I tackle some pressing deadlines.

miaoyinb avatar Apr 23 '25 14:04 miaoyinb

Added a unit test. It should work with the this PR's fix but will fail if run with the current devel

Not sure why some civet tests failed (not familiar with libmesh's testing system).

miaoyinb avatar Apr 28 '25 17:04 miaoyinb

Your CPPUNIT_TEST() invocation is in the #ifdef LIBMESH_HAVE_POLY2TRI block, but you defined it in the middle of the #ifdef LIBMESH_HAVE_TRIANGLE block. For builds that have Poly2Tri enabled but not Triangle (which is common and includes all MOOSE builds - Triangle is GPL-only), the unit tests want to include your new test but don't see a definition for it.

The ideal thing to do for tests that can be supported by both Poly2Tri and Triangle backends is to create a single implementation that works with a reference-to-base-class and two instantiations that each pass in one of the different subclasses. There's plenty of examples in there to look at.

roystgnr avatar Apr 28 '25 21:04 roystgnr

I added the unit test as directed. thanks.

miaoyinb avatar Apr 30 '25 03:04 miaoyinb

@roystgnr just a friendly reminder, thanks.

miaoyinb avatar Jun 04 '25 02:06 miaoyinb

Hi @roystgnr , can you take a look at this again. The unit test has been added so that you can check when the node ids are changed. Thanks.

miaoyinb avatar Jun 10 '25 15:06 miaoyinb