Polygon implementation for modern (spherical mean value) interpolation.

[lewisn-met]

Description

This PR relates to issue ecmwf/atlas#328 and adds:

• A SphericalPolygon3D class as an alternative interpolation method to the Triag3D and Quad3D classes currently used in the finite-element interpolation method (FiniteElement.h). This computes weights using the spherical mean value formula rather than the bilinear/barycentric formulae currently implemented.

• A test suite for the new class, in test_SphericalPolygon3D.cc, which validates that the interpolation method correctly computes the areas, normals, and weights of an example triangle, planar quadrilateral, and non-planar quadrilateral on the surface of the unit sphere.

(Draft) PR integrating this class into the existing set of interpolation methods in atlas to follow, with the exact implementation to be fixed after further discussion with myself and @odlomax.

Contributor Declaration

By opening this pull request, I affirm the following:

• All authors agree to the Contributor License Agreement.
• The code follows the project's coding standards.
• I have performed self-review and added comments where needed.
• I have added or updated tests to verify that my changes are effective and functional.
• I have run all existing tests and confirmed they pass.

[wdeconinck]

Thanks @lewisn-met for this contribution. I am wondering in how far there is overlap with the ConvexSphericalPolygon which is used in the "ConservativeSphericalPolygon" interpolation, and if missing functionality can be merged perhaps. @sbrdar could you also have a look please?

[codecov-commenter]

Codecov Report

:x: Patch coverage is 84.61538% with 32 lines in your changes missing coverage. Please review. :white_check_mark: Project coverage is 79.78%. Comparing base (4a5d4e9) to head (70ce3ff). :warning: Report is 2 commits behind head on develop.

Files with missing lines	Patch %	Lines
src/tests/interpolation/test_SphericalPolygon3D.cc	80.00%	32 Missing :warning:

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #329      +/-   ##
===========================================
+ Coverage    79.27%   79.78%   +0.50%     
===========================================
  Files          909      909              
  Lines        61597    70950    +9353     
===========================================
+ Hits         48832    56605    +7773     
- Misses       12765    14345    +1580     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[lewisn-met]

Thanks for taking a look at this and for the kind words @wdeconinck. From a brief read of the ConvexSphericalPolygon class, the method used to check for the left of a great circle is precisely what I used to check if a point was inside SphericalPolygon3D. As for the intersection method, I will need to look a bit further, but they seem to be distinct.

[lewisn-met]

On a slightly unrelated note, I notice that the linux intel-classic ci and gnu-14 ci are are currently failing for this PR. These failures seem to be related to the pluto-benchmark and either the SphericalVector or Cubic2D classes. Is this something you'd like me to take a look at?

[wdeconinck]

On a slightly unrelated note, I notice that the linux intel-classic ci and gnu-14 ci are are currently failing for this PR. These failures seem to be related to the pluto-benchmark and either the SphericalVector or Cubic2D classes. Is this something you'd like me to take a look at?

My hunch is that it is due to an updated Eigen version 5.0!!! that happened under the hood without us realising: linux intel-classic:

  -- Found Eigen3: /home/linuxbrew/.linuxbrew/share/eigen3/cmake/Eigen3Config.cmake (found version "5.0.0")

Any help to fix this is always welcome of course.

[wdeconinck]

I have followed up on my hunch, and reverting the CI to use "eigen@3" seems to fix the compilation. I have added this fix to develop and a 0.44.1 hotfix tag.

This is not a proper fix for compatibility Eigen 5.0.0 but perhaps these are infancy issues that will be sorted upstream in Eigen in due course. At least we can work without failing CI now. Please rebase on develop.

[lewisn-met]

I have now merged the fix from develop into this PR. Thanks for following your hunch so promptly - fingers crossed that any issues with Eigen 5. get sorted out soon!

I also noticed that the validate method included in your ConvexSphericalPolygon class nicely verifies the polygon is on the sphere and points are not too close. I could not see anywhere in the FiniteElement class that does a similar check, so it may be worth sharing this functionality between the methods.

[sbrdar]

Hi @lewisn-met. Many thanks for this interesting development. It seems to me that SphericalPolygon3D can be class-derived from ConvexSphericalPolygon. Can we assume that SphericalPolygon3D has to be convex spherical polygon?

[lewisn-met]

Hi @sbrdar, thank you for taking a look at this. Indeed, we can assume SphericalPolygon3D is convex, though it need not be (the weights are still positive for star-shaped polygons and well-defined but negative otherwise).

With regards to merging functionality, incorporating the methods of SphericalPolygon3D into ConvexSphericalPolygon, either directly or via inheritance, is certainly an option. The main differences/obstacles here are the representation of the polygon (SphericalPolygon3D requires the triple products from your inLeftHemisphere() function to be stored in addition to the vertices) and encapsulating what is actually necessary for each use case.

Hopefully we can discuss this soon.

[sbrdar]

Hi @lewisn-met I am thinking that the Floater's barycentre could be used for the 2nd order cons. remapping in Atlas, so I cherry-picked your implementation with the unit test on SphericalVector3D, and added it to ConvexSphericalPolygon, see the usage here. Maybe something to discuss soon.

[lewisn-met]

Would incorporating this method into the cons. remapping be in addition or as an alternative to using Floater's interpolation in the finite-element class?

Our original motivation for the algorithm was as a drop-in alternative to the Quad3D and Triag3D classes used for "continuous" interpolation, see the implementation here. However, I can see how Floater's ideas could be useful in both settings and I don't expect the use of ConvexSphericalPolygon in FiniteElement.cc to be too different from the branch I have linked.