meshio
meshio copied to clipboard
nschloe
OpenFOAM file format
Hi,
I'm trying to implement support for exporting meshes into the OpenFOAM polyMesh format. This format consists of a directory with several files describing the mesh. This would permit bypassing the outdated gmsh 2.2 format, as OpenFOAM does not support later gmsh formats.
For now I've added basic write support, but physical groups are still missing. These will be required for most use cases as OpenFOAM expects outer mesh surfaces to be named so that appropriate boundary conditions can be assigned. However, I'm struggling to understand how these are defined in Mesh objects (specifically for those exported by pygmsh). Any help on this would be appreciated.
I'll keep working on this. Creating pull request now to see if there is any interest/suggestions.
This needs tests. Take a look at the other format tests -- they essentially all so the same.
I've added named patches, which should mean that writing functionality is complete for basic use cases. I'll start writing tests.
Is the ability to read OpenFOAM meshes required for merging?
This looks good. How about adding a test file? Doesn't need to be fancy. The most important thing is the write-read test; take a look at the other tests, they all look more or less the same.
@nschloe Yeah absolutely, I'm just new to writing tests so I need to do some reading.
I'd check out the other tests, some are only 5 lines, e.g., dolfin or tetgen. You don't need more than that.
ok, I've added simple tests but they fail because the reader doesn't always output the correct node ordering. If I make the check ignore node order the tests pass, so I think that's the only issue.
Will need to add code to reorder the nodes. Not sure how difficult that will be. OpenFOAM defines elements by their faces in any order so this may take some thinking to get right.
The error is
ModuleNotFoundError: No module named 'helpers'
You'll also need to run make format to make the linter happy.
My bad! I was running pytest from inside the tests folder (couldn't use tox because it complains that numpy is missing for some reason?) so I wasn't getting that error. Should be fixed now.
Edit: If anyone else is having this issue running tox -r just fixed it
ok, I've added simple tests but they fail because the reader doesn't always output the correct node ordering.
Right, that's the issue now. Would be cool to get this merged!
hmm ok there's a complication: I don't think there's any way of ensuring the same node order in a write-read cycle because the 'rotation' of elements is lost in the way OpenFOAM stores mesh data. I could make it respect the node ordering as given in the wiki but I can't use the helpers.write_read function in its current state because it is too restrictive.
Seems to me the options are:
- Change the
helpers.write_readfunction to allow any equivalent but non-identical node ordering (sounds very difficult to me) - Write a more basic test which doesn't check node ordering (easy but then the test isn't too useful)
Any thoughts?
'rotation' of elements is lost in the way OpenFOAM stores mesh data.
Well there must be some ordering, right? The quad [0,1,2,3] is different from [0,1,3,2]. How does openfoam store those?
Well I think the issue is that the quad [0, 1, 2, 3] is in a way the same as [3, 0, 1, 2].
OpenFOAM only deals with 3D elements, but it stores them as a list of (right-hand rule) faces in any order. So for example if we take a hexahedron, it is possible to call any vertex 'vertex 0' and then number the rest of the points in a way which respects the wiki convention. I don't think there's any way of ensuring that the same vertex is selected as 'vertex 0'.
Tests work now. However I had to duplicate a lot of code from helpers.write_read to make a version that's less restrictive. Might be better to add the changes as an optional argument to the original function instead?
I still don't like very much the fact that you don't get the same input/output from a meshio cycle. One suggestion would be to normalize the faces when reading them, i.e., from the hex
[0, 1, 2, 3, 4, 5, 6, 7]
you write
[
[3, 2, 1, 0],
[0, 1, 5, 4],
[1, 2, 6, 5],
[2, 3, 7, 6],
[3, 0, 4, 7],
[4, 5, 6, 7],
]
When reading, you normalize the faces data in the following way:
- You take the first row, revert.
- In the rest of the faces, you look for the face which has no points in common with the first row.
- Then you look for the two other rows which have have the first index of the reverted first row (
0). They have one other index in common (4) - Roll the row you found in step 2 such that that index comes first.
- Concatenate the two rows.
You end up with
[0, 1, 2, 3, 4, 5, 6, 7]
This should read all meshes correctly, and return the same mesh from a meshio write-read cycle.
What do you think?
This is exactly the method I was trying to implement... except I hadn't realized the need to ensure the second face starts at the right index! I'll fix that.
My concern was that I think there are still variations possible because the index assigned to each point is not necessarily preserved between read and write. So for example if we rotate a hex by 90deg along any axis, the node order is still good but the point that was previously 0 might now be called 4 etc.
Having said that, I've just realized that I was unnecessarily rearranging the order of faces on write because I misread the docs (I didn't see that I could skip points in the neighbour file using index -1, so I was sorting all neighboured faces together 🤦). If the faces were kept in the same order it should be possible to preserve point indices?
My concern was that I think there are still variations possible because the index assigned to each point is not necessarily preserved between read and write. So for example if we rotate a hex by 90deg along any axis, the node order is still good but the point that was previously 0 might now be called 4 etc.
How would that happen if you start from the exact same input data? If the order of the input mesh indices is different, it's okay if the output is also in a different order. The important thing is that input->write->read returns the input mesh.
I plan to make a new release soon. Are you still interested in getting this to work?
Yeah I am just have less time now because I started a job. I'll get around to it.
I'm interested in seeing this PR merged and I'm willing to help if @nicola-sorace is too busy right now. What needs to be done still?
@NauticalMile64 Read and write should work at least at a basic level, but we ran into an issue writing tests because there doesn't seem to be a way of providing an arbitrary node order in the OpenFOAM mesh format. Don't quite remember the specifics anymore but I think the main obstacle was something like elements which are part of a named group have to be sequential. As a result, the standard test case of 'write->read->mesh is identical' doesn't work because the element order comes back different.
Solution is either to find some way of preserving node order in the write->read process, or to write a new test case which ignores element order. The 'specialty function' mentioned above is doing the latter, but I copy-pasted it from the standard test and changed it slightly to ignore order; this is a lot of code duplication which is not great.
Personally I think the simplest/cleanest solution would be to add a parameter to the original test case for ignoring element order (and maybe add a warning for file formats which don't preserve order?). There might be better solutions though.
@nschloe @NauticalMile64 Hello,
Decided to try this again. Surprisingly, I removed the specialty function and all the tests are passing now. Maybe the write_read() function has been "fixed" since last time, for some other file format?
oh never mind, I updated the code wrong and it was just skipping all cells. Node order issue still exists.
Less the tests issue, does this in theory convert the mesh to other formats retaining boundaries and regions/blocks labels?
Unless something has broken in recent releases, yes. There are probably some limitations though, for example I think it only supports labels for faces and not for cells. And I think a few cell types haven't been implemented yet. But for my use-case it was working.