Edges position computed by `MpasMeshConverter.x` are incorrectly not placed at primal and dual meshes intersections

[favba]

I'm using the MpasMeshConverter.x tool to create planar periodic meshes by feeding cells and vertices position, together with the connectivity, and I noticed that the area of the dual cells (triangles) computed by the tool is not correct.

I'll provide a MWE when I have the time. Meanwhile, this is the ncdump of the mesh created where I noticed this issue mesh_distorted_dump.txt The mesh was created based on the xCell,yCell,zCell,xVertex,yVertex,zVertex, cellsOnVertex and meshDensity info only. The rest of the data was computed by MpasMeshConverter.x. I haven't checked yet if the other areas are correct.

[xylar]

@favba, this is much appreciated. It may relate to https://github.com/MPAS-Dev/MPAS-Tools/issues/289.

[xylar]

Just by inspection, I don't see anything wrong with the computation of areas so this will need some deeper debugging.

[favba]

An important observation: The meshes I'm generating are purposely inhomogeneous. I first generate a homogeneous mesh with planar_hex, then I disturb the voronoi cells generator points so I don't have a perfectly homogeneous mesh anymore, then recompute the triangles circumcenters based on the new Voronoi Cells position. Maybe the tool is assuming a homogeneous mesh?

I also noticed that 1) while each triangle area is incorrect, the sum of all triangle areas still gives the correct result. 2) The Voronoi Cells areas are correct 3) The kiteAreasOnVertex are also incorrect (but consistent with the incorrect triangle area)

Here is a small mesh input where I see the issue:

netcdf test_mesh_distorted_incomplete {
dimensions:
    nCells = 16 ;
    nEdges = 48 ;
    nVertices = 32 ;
    maxEdges = 6 ;
    TWO = 2 ;
    maxEdges2 = 12 ;
    vertexDegree = 3 ;
variables:
    double xCell(nCells) ;
    double yCell(nCells) ;
    double zCell(nCells) ;
    double xVertex(nVertices) ;
    double yVertex(nVertices) ;
    double zVertex(nVertices) ;
    int cellsOnVertex(nVertices, vertexDegree) ;
    double meshDensity(nCells) ;

// global attributes:
        :is_periodic = "YES" ;
        :x_period = 4. ;
        :y_period = 3.46410161513775 ;
        :dc = 1. ;
        :nx = 4 ;
        :ny = 4 ;
        :on_a_sphere = "NO" ;
        :sphere_radius = 0. ;
        :history = "Fri Apr 26 09:28:30 2024: /home/user/.julia/environments/mpas/.CondaPkg/env/bin/planar_hex --nx 4 --ny 4 --dc 1.0 -o test_mesh.nc" ;
data:

 xCell = 0.419612, 1.558266, 2.58442, 3.58442, 0.904384, 2.03996, 3.01458, 
    3.947584, 0.554684, 1.47876, 2.52916, 3.559976, 1.06237, 1.955774, 
    3.099324, 4.031104 ;

 yCell = 0.978931247105323, 1.0035951433495, 0.959562902147395, 
    0.952151030632665, 1.84684601889757, 1.85911550919761, 1.89339337355814, 
    1.88723780104449, 2.7657861399324, 2.70213053562675, 2.69079910173104, 
    2.70477503713879, 3.63103041163494, 3.62181813590841, 3.51498680968391, 
    3.55619060027713 ;

 zCell = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;

 xVertex = -0.00975376456278611, 0.0128193185794142, 0.983698434302909, 
    0.993667963873752, 2.08376564654827, 2.06025408134221, 3.08692684885916, 
    3.08208886589458, 0.439874804893943, 0.412690218024825, 1.46963538486421, 
    1.47475428549263, 2.5180783205772, 2.53886634136681, 3.48291129552627, 
    3.47871316725034, 0.0414356733004868, 0.0749339411046957, 
    1.0393849077105, 0.994201634447349, 2.00735865125447, 2.00121085131072, 
    3.0412042843806, 3.04779109408958, 0.526858412042012, 0.565049572904608, 
    1.51219928594855, 1.50541588331441, 2.54952770613702, 2.50410416278026, 
    3.5511795753388, 3.57780005166295 ;

 yVertex = 1.33260363870946, 0.628619095681241, 1.23320352104529, 
    0.772941905785019, 1.27108389686055, 0.723156150650342, 1.2940777321148, 
    0.64134340516813, 2.19608711214069, 1.55213897878609, 2.08775166316379, 
    1.61398294434437, 2.13760078737025, 1.54653651892945, 2.1675830615144, 
    1.53127011657355, 2.99441725708763, 2.44827084364361, 3.0629513515751, 
    2.40703458070538, 3.01151261497334, 2.44162453899308, 2.94588294128612, 
    2.46006294864359, 3.86753038251839, 3.34126976021641, 3.92970761640241, 
    3.27185506977303, 3.80366819766577, 3.31744290363107, 3.85296232576967, 
    3.25096850257545 ;

 zVertex = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;

 cellsOnVertex =
  1, 8, 4,
  4, 16, 1,
  2, 5, 1,
  1, 13, 2,
  3, 6, 2,
  2, 14, 3,
  4, 7, 3,
  3, 15, 4,
  5, 9, 8,
  8, 1, 5,
  6, 10, 5,
  5, 2, 6,
  7, 11, 6,
  6, 3, 7,
  8, 12, 7,
  7, 4, 8,
  9, 16, 12,
  12, 8, 9,
  10, 13, 9,
  9, 5, 10,
  11, 14, 10,
  10, 6, 11,
  12, 15, 11,
  11, 7, 12,
  13, 1, 16,
  16, 9, 13,
  14, 2, 13,
  13, 10, 14,
  15, 3, 14,
  14, 11, 15,
  16, 4, 15,
  15, 12, 16 ;

 meshDensity = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ;
}

This is what it looks like:

After feeding this to MpasMeshConverter.x I get the following test_mesh_distorted.txt

Take the 13th vertex in the output mesh for example. The cells (Fortran) indices surrounding it are (7, 11, 6). The position of such cells are:

((3.01458, 1.8933933735581412),
 (2.52916, 2.6907991017310433), 
 (2.0399599999999998, 1.8591155091976144))

A triangle formed by such vertices has an area of ~0.3969033658548805, but the area computed by MpasMeshConverter.x was 0.43357495536806157

[xylar]

Thanks @favba. Can you also attach the NetCDF file itself (before conversation)?

[favba]

Here is a zip file with the three meshes: test_meshes.zip

test_mesh.nc is the unaltered mesh produced by planar_hex test_mesh_distorted_incomplete.nc is the incomplete distorted mesh which is constructed based on test_mesh.nc and then fed to MpasMeshConverter.x test_mesh_distorted_hdf5.nc is the final one, with the wrong triangle areas.

[xylar]

@favba, where are you getting the MpasMeshConverter.x from? From the conda package? Or building it yourself?

[favba]

From the conda package:

 mpas_tools v0.33.0 (channel=conda-forge)

[xylar]

@favba, also, have you verified that xEdge and yEdge are what you think they should be? It seems to me like the most plausible cause of the area mismatch is that xEdge and yEdge are not in the center of the line connecting cell centers, which might indicate that the perturbed mesh isn't quite a Voronoi mesh but I haven't checked this.

[xylar]

Thanks, that's helpful. There are 2 versions of the mesh converter code and that tells me which one to look at.

[xylar]

@favba, as I suspected, it looks like the perturbed mesh you're supplying isn't quite a valid Voronoi mesh The edge centers are at the centers of lines connecting vertices but they're not at the centers of lines connecting the cell centers:

This is for the 13th vertex. Edge centers are the red dots, cell centers in blue and the vertex in black. You can see that the edge centers aren't on the lines as they should be. The mesh converter doesn't enforce this property but it very much assumes it to be true.

[favba]

Thanks, @xylar, I did not check that! I was indeed expecting the behavior you were describing. Note that I don't inform the edge positions in the mesh I feed to MpasMeshConverter.x and hoped it computed the one expected by MPAS to me...

I believe the mesh I am supplying is still a voronoi mesh: any point in the (primal cell) edge should be equidistant to the cell points generators. It just happen that my generators points are (purposely) not in optimized positions. But I'll double-check on that again.

One of the reasons I'm building this kind of distorted mesh is precisely to investigate the effects of this issue: not having {x,y,z}Edge fall in the edge midpoint. This issue does occur on real variable resolution MPAS meshes... I was hoping to investigate the effects of that using simpler 2D planar cases.

It seems to me this is a bug with the mesh converter then, is it not? It should compute {x,y,z}Edge as the midpoint of the dual mesh edge, not of the primal mesh edge. And, if it is indeed a Voronoi mesh, that would fall in the primal and dual meshes intersection (but not necessarily in the primal cell edge midpoint).

[xylar]

We would need to discuss with a broader team how we expect the converter to behave in cases like this. My understanding was that this was the expected behavior -- edge locations are computed to be halfway between vertices and expected to be exactly half way between cell centers. Thus, I think the converter is doing what it's designed to do. I would need to really scrutinize the mesh spec https://mpas-dev.github.io/files/documents/MPAS-MeshSpec.pdf to see if it should be the other way around but we've been using the current approach for about 15 years now...

[xylar]

In Fig. 3.1, the edge locations are depicted as being halfway between vertices.

[xylar]

One of the reasons I'm building this kind of distorted mesh is precisely to investigate the effects of this issue: not having {x,y,z}Edge fall in the edge midpoint. This issue does occur on real variable resolution MPAS meshes... I was hoping to investigate the effects of that using simpler 2D planar cases.

My understanding is that, at least the variable resolution meshes we are using from the JIGSAW tool, are true Voronoi grids where the centers of edges are both halfway between cell centers and vertices. I know our group (working mostly on MPAS-Ocean) had considered alternatives where that was not the case but I don't think we had revisited the mesh converter in that process. Presumably, for meshes that aren't strictly Voronoi, one would need to have a different tool for computing areas depending on your needs. There would be 2 sets of kite areas depending on where you define edge centers.

[matthewhoffman]

@favba , if you have manually modified the cell center locations but kept the triangulation topology unchanged, you could use this tool before running MpasMeshConverter.x to generate a valid Voronoi mesh from your pointset: https://github.com/MPAS-Dev/MPAS-Tools/tree/master/mesh_tools/points-mpas I believe that you will get the correct behavior from MpasMeshConverter.x if you've first processed your pointset with that intermediate tool.

[mwarusz]

My understanding is that, at least the variable resolution meshes we are using from the JIGSAW tool, are true Voronoi grids where the centers of edges are both halfway between cell centers and vertices.

@xylar I don't think that's correct. This came up when I was doing convergence testing for Omega. The only property that a Voronoi grid guarantees is that the edge points are halfway between cell centers. My understanding is that what you call a true Voronoi grid is only possible for uniform planar meshes. In particular, edge points on spherical meshes generated by Jigsaw are not halfway between vertices.

I believe this issue is a genuine bug in MPAS mesh converter. The figure you created was very helpful in tracking it down. As you can see here: https://github.com/MPAS-Dev/MPAS-Tools/blob/d29c4f218dbd9a9ec4b3bc5de8e797fb1c0e1d8a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp#L1253 and here: https://github.com/MPAS-Dev/MPAS-Tools/blob/d29c4f218dbd9a9ec4b3bc5de8e797fb1c0e1d8a/mesh_tools/mesh_conversion_tools_netcdf_c/pnt.h#L374-L378

the converter intends to put edge points on the intersection of the line connecting the centers with the line connecting the vertices. So it doesn't make sense for such point to not lie on a triangle edge, if the intersection is computed correctly.

This lead me to suspect that the existing formula for the intersection is not correct. In #567 I simply replaced the planar line-line intersection formula with one taken from Wikipedia. Running this modified mesh converter on the mesh that @favba supplied produces the correct area for the triangle corresponding to the 13th vertex.

[favba]

Thanks, @mwarusz ! I agree with everything stated. I've updated the issue title to reflect what seems to be the source of the problem. Also, for the fix, if we expect the input to be a valid Voronoi diagram (which is a necessary condition for MPAS meshes), then the Edge position in the plane can be computed using a simpler formula: 0.5*(cell1 + cell2) , where cell1 and cell2 are the coordinates vectors of the Voronoi cells separated by the given edge.

[favba]

Thanks for the suggestion, @matthewhoffman . But, reading the Readme of that utility, it seems it only works for meshes on the sphere, no?

[xylar]

Thanks @favba for persisting on this and @mwarusz for weighing in. Since this tool is fundamental to the definition of MPAS meshes, we won't be able to make things change until we carefully consult with all he stakeholders.

[xylar]

Addressed by #567. Thank you @favba and @mwarusz!