From 293b8a4a00c54af77dbc0c7f18122318000a12d8 Mon Sep 17 00:00:00 2001 From: Xylar Asay-Davis Date: Thu, 24 Sep 2026 07:32:08 -0500 Subject: [PATCH] Fix law-of-cosines typo in triangleAngleQuality The second angle of each dual triangle in buildMeshQualities() used b_len * c_len where the law of cosines needs b_len * b_len, so triangleAngleQuality (and potentially obtuseTriangle) were wrong wherever the triangle's b and c sides differ. Fix both the netcdf_c converter (built into the conda package) and the legacy converter. Add a test that recomputes the dual-triangle angles from dcEdge and checks triangleAngleQuality, obtuseTriangle and that the angles sum to pi. Co-Authored-By: Claude Opus 5.5 (1M context) --- conda_package/tests/test_conversion.py | 36 +++++++++++++++++++ .../mpas_mesh_converter.cpp | 2 +- .../mpas_mesh_converter.cpp | 2 +- 3 files changed, 38 insertions(+), 2 deletions(-) diff --git a/conda_package/tests/test_conversion.py b/conda_package/tests/test_conversion.py index 3a02eb98e..23377d622 100755 --- a/conda_package/tests/test_conversion.py +++ b/conda_package/tests/test_conversion.py @@ -47,6 +47,41 @@ def test_conversion_angle_edge(): assert np.max(np.abs(angle_diff)) < 1.0e-10 +def test_conversion_triangle_angle_quality(): + ds_mesh = xarray.open_dataset( + get_test_data_file('mesh.QU.1920km.151026.nc') + ) + ds_mesh = convert(dsIn=ds_mesh) + + edges_on_vertex = ds_mesh.edgesOnVertex.values - 1 + assert np.all(edges_on_vertex >= 0) + dc_edge = ds_mesh.dcEdge.values + a_len = dc_edge[edges_on_vertex[:, 0]] + b_len = dc_edge[edges_on_vertex[:, 1]] + c_len = dc_edge[edges_on_vertex[:, 2]] + + # law of cosines for the angle opposite each side of the dual triangle + angle1 = np.arccos( + np.clip((b_len**2 + c_len**2 - a_len**2) / (2 * b_len * c_len), -1, 1) + ) + angle2 = np.arccos( + np.clip((a_len**2 + c_len**2 - b_len**2) / (2 * a_len * c_len), -1, 1) + ) + angle3 = np.arccos( + np.clip((a_len**2 + b_len**2 - c_len**2) / (2 * a_len * b_len), -1, 1) + ) + angles = np.stack([angle1, angle2, angle3], axis=1) + assert np.max(np.abs(angles.sum(axis=1) - np.pi)) < 1.0e-10 + + min_angle = angles.min(axis=1) + max_angle = angles.max(axis=1) + quality = ds_mesh.triangleAngleQuality.values + assert np.max(np.abs(quality - min_angle / max_angle)) < 1.0e-10 + + obtuse = (max_angle > 0.5 * np.pi).astype(int) + assert np.array_equal(ds_mesh.obtuseTriangle.values, obtuse) + + def test_masks_to_int_dataset_copy(): ds_in = xarray.Dataset( data_vars={ @@ -71,3 +106,4 @@ def test_masks_to_int_dataset_copy(): if __name__ == '__main__': test_conversion() test_conversion_angle_edge() + test_conversion_triangle_angle_quality() diff --git a/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp b/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp index 83d052aea..a578eb0b0 100644 --- a/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp +++ b/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp @@ -2314,7 +2314,7 @@ int buildMeshQualities(){/*{{{*/ #endif angle1 = acos( max(-1.0, min(1.0, (b_len * b_len + c_len * c_len - a_len * a_len) / (2 * b_len * c_len)))); - angle2 = acos( max(-1.0, min(1.0, (a_len * a_len + c_len * c_len - b_len * c_len) / (2 * a_len * c_len)))); + angle2 = acos( max(-1.0, min(1.0, (a_len * a_len + c_len * c_len - b_len * b_len) / (2 * a_len * c_len)))); angle3 = acos( max(-1.0, min(1.0, (a_len * a_len + b_len * b_len - c_len * c_len) / (2 * a_len * b_len)))); minAngle = min(angle1, min(angle2, angle3)); diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp index 6f38e7424..aafeda8c4 100755 --- a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp +++ b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp @@ -2395,7 +2395,7 @@ int buildMeshQualities(){/*{{{*/ angle1 = acos( max(-1.0, min(1.0, (b_len * b_len + c_len * c_len - a_len * a_len) / (2 * b_len * c_len)))); angle2 = acos( max(-1.0, min(1.0, - (a_len * a_len + c_len * c_len - b_len * c_len) / (2 * a_len * c_len)))); + (a_len * a_len + c_len * c_len - b_len * b_len) / (2 * a_len * c_len)))); angle3 = acos( max(-1.0, min(1.0, (a_len * a_len + b_len * b_len - c_len * c_len) / (2 * a_len * b_len))));