Mesh geometry utils

[pabloyesteb]

Adding functions to compute the cells connectivity and edge normals of a mesh. This is the same functionality that was intended to be added with the issue #68 which was never closed and can now be discarded. The reason is that new algorithms have been tried which reduce the computation time of the cells connectivity from 4 hours to around 20 seconds for a mesh with ~1e6 elements.

The basic idea of the algorithm is to create a dictionary of edges. For every edge, the dictionary stores the (two at most) cells that contain that given edge. Computing the cells connectivity is then trivial and computationally inexpensive (I believe the algorithm is o(n)).

The two functions being implemented are:
compute_cells_connectivity(conec)
returns dict associating a list of neighbors to every cell index

compute_edge_normals(xyz, conec)
    returns dict associating a list of edge normals to every cell index.

conec should be an array containing the nodes defining each cell in the mesh (this is called self.connectivity in the Mesh class)
xyz should be an array containing the coordinates of the nodes (this is called self.xyz in the Mesh class)

The two functions are to be implemented in pyLOM/vmmath/wrapper.py

P.S.
The algorithm for computing the cells connectivity can easily be used to ensure coherence in the direction of surface normals for open meshes (all of them pointing either outwards or inwards). I believe the actual function for calculating the normals in vmmath.wrapper does not ensure coherence in the direction of the normals. This can be modified in a future issue under request.

[ArnauMiro]

@pabloyesteb can you create a branch linked to this issue and push your changes? The normal computation code relies on the correct orientation of the elements so I would say is coherent. I'm open to explore how to ensure coherence with your new algorithm, that could be done in this same issue.

I will then close #68 and delete its associated branch.

[ArnauMiro]

@pabloyesteb I'll review your PR in the following days and approve it. Looks good on first glance, just need to make sure that we are consistent with the compiled version. The example uses any of the CETACEO meshes? Could we, for example, provide a test for the cylinder case, as it is more straightforward to the users?

[pabloyesteb]

Thanks. The example should work for any (2d) mesh file that can be read by the mesh class. As of now, the user needs to provide such file to run the example. I'll check the cylinder case and try to include it in the example.

[ArnauMiro]

Can you check that it works for the cylinder example? I will use that to adapt your code to standard C.

[pabloyesteb]

Sure. Could you please provide me with the cylinder h5 file? I can't find it anywhere.

[ArnauMiro]

You can get them from here: https://huggingface.co/datasets/bef-18/pyLOM_examples/tree/main

[pabloyesteb]

I think I got it right. The code works for sure, and looking at the figures, it all seems to be right. I wasn't able to visualize the results in ParaView though. I tried building a vtu file from the original h5 file that ParaView can read but failed. Do you have an vtkhdf writer? How did you get the vtk.hdf files from cetaceo?

[ArnauMiro]

@pabloyesteb no need to go for the VTU, we have VTKHDF already implemented. Check the writer in inp_out/io_paraview.py in particular:

pv_writer(Mesh,Dataset,casestr,basedir='./',idim=0,instants=[0],times=[0.],vars=[],fmt='vtkh5',mode='w')

[pabloyesteb]

Alright. I'm adding new functionalities to the geometry routines I need for the gnn implementation and when I'm done I'll check them with the pv_writer.

[ArnauMiro]

Also do a pull, I've adapted some of your code to the compiled version.