Skip to content

fix the quadratic element issue in TriangulatorInterface #4145

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: devel
Choose a base branch
from

Conversation

miaoyinb
Copy link
Contributor

closes #4144

@miaoyinb
Copy link
Contributor Author

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

@moosebuild
Copy link

moosebuild commented Apr 18, 2025

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.

@miaoyinb miaoyinb marked this pull request as ready for review April 19, 2025 03:24
Copy link
Member

@roystgnr roystgnr left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@miaoyinb
Copy link
Contributor Author

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.

@roystgnr
Copy link
Member

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
Copy link
Member

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.

@miaoyinb
Copy link
Contributor Author

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

@miaoyinb miaoyinb force-pushed the quadraticTriangulatorFix branch 4 times, most recently from 6922c7f to 238f4cd Compare April 28, 2025 16:01
@miaoyinb
Copy link
Contributor Author

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).

@roystgnr
Copy link
Member

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.

@miaoyinb miaoyinb force-pushed the quadraticTriangulatorFix branch 2 times, most recently from b7d819c to 216ba86 Compare April 29, 2025 14:23
@miaoyinb miaoyinb force-pushed the quadraticTriangulatorFix branch from 216ba86 to 33866f9 Compare April 29, 2025 14:30
@miaoyinb
Copy link
Contributor Author

I added the unit test as directed. thanks.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Incorrect Outer Boundary Midpoints in TriangulatorInterface::increase_triangle_order()
3 participants