nschloe / meshio

:spider_web: input/output for many mesh formats
MIT License
1.93k stars 399 forks source link

MED 4.0 to .msh or .xml with entities #347

Closed RemDelaporteMathurin closed 5 years ago

RemDelaporteMathurin commented 5 years ago

Hi all,

I wish to convert a mesh from SALOME (MED 4.0) to .msh or .xml alongside with cell attributes.

Here's an exemple of a mesh from SALOME with 1 volume marked. MED.zip

Meshio seems to read it without issues. But, when converting to .msh and after opening it in GMSH, here's the mesh file I get. gmsh.zip

No sign of the entity "tungsten".

Thanks guys for your help!

Rem

PS: Screenshot of the med file re-imported in SALOME, the group is there as expected. image

RemDelaporteMathurin commented 5 years ago

Update: I tried to convert the MED file to another MED file with meshio, it seems it removed the topological information ... I'm using SALOME 9.2.2 .

image

gdmcbain commented 5 years ago

Earlier this month we started putting together a test suite of authoritatively generated mesh files #335, e.g. MSH by Gmsh #342. Could this one serve as such an exemplar for MED? I don't know much about MED myself; is Salome its authoritative generator?

RemDelaporteMathurin commented 5 years ago

I think it is though I'm far from being an expert. Anyway, this one can definitely serve a an exemple.

RemDelaporteMathurin commented 5 years ago

I noticed in med_io.py, the function read(filename) doesn't seem to take into account groups which should be under FAS. Moreover in write(filename, mesh, add_global_ids=True), the Subgroups seem not to be filled (just the blank fields). The groups seem not be read nor written in any case. This may be the reason why the groups vanish when converting from MED to MED ?

@nschloe Please correct me if I'm wrong I might have misunderstood something there

tianyikillua commented 5 years ago

@RemiTheWarrior That's correct! I'm not sure if there are already established conventions as how to store node/element groups in a meshio-mesh class...Using these information to write to a MED file should be easy...

gdmcbain commented 5 years ago

Thank you. I won't be able to do too much with this myself over the weekend, away from the office, but I wonder whether this is associated with the conversion of tagged subsets of elements #290 between formats. A neutral representation of these is yet to be established.

RemDelaporteMathurin commented 5 years ago

@tianyikillua I managed to find some info on MED format here but that's not really precise.

tianyikillua commented 5 years ago

The best way to learn MED is to open MED files using a HDF5 viewer, and happy hacking...I used this way to implement read_data and write_data for MED files...

Concerning MED to MED while reading and writing available node/element groups, it should be possible to do so, in the first time, within the med module, as

A neutral representation of these is yet to be established.

If you are interested feel free. Otherwise I may find some time next week...

RemDelaporteMathurin commented 5 years ago

Currently downloading HDFView 😉 This is kinda new to me, but I'll definitely try to look into it

RemDelaporteMathurin commented 5 years ago

Hi all,

Just to say that I found a workaround using for converting .med to .xml using MED 3.2: 1) Convert from .med to .msh by opening it in GMSH 2) Convert from .msh to .xml with dolfin-convert

Note that 2) creates as expected 3 .xml files (mesh.xml alongside mesh_facet_region.xml and mesh_physical_region.xml) whereas meshio-convert doesn't and returns this warning: WARNING:root:DOLFIN XML can only handle one cell type at a time. Using tetra, discarding triangle.

Here are the files. MED 3.2 generated with SALOME 8.3.0

.msh converted with GMSH 3.0.0

.xml converted with dolfin-convert

gdmcbain commented 5 years ago

Thanks @RemiTheWarrior, this is really useful; I propose to include these into the test suite so that we can test that meshio can correctly perform these conversions too.

gdmcbain commented 5 years ago

@RemiTheWarrior, I notice that Mesh_1_830.med contains 72 2-node line elements which are not in Mesh_1_830.msh. Are they significant? Both contain 270 3-coordinate points (matching to within ±1/10¹³), 438 3-node triangles (identical connectivity), and 839 4-node tetrahedra (same connectivity but with middle two points permuted).

The main difference is of course that the meshio.Mesh read from the MED contains no cell_data or field_data whereas that from your MSH has field_data:

{'top': array([1, 2]), 'bottom': array([2, 2]), 'walls': array([3, 2]), 'tungsten': array([4, 3])}

and cell_data[t]['gmsh:physical'] for t = 'triangle' (438 values in {1, 2, 3}) and tetra (839 values, all 4).

But I just wanted to check about the line-elements before proceeding.

gdmcbain commented 5 years ago

The best way to learn MED is to open MED files using a HDF5 viewer

Another good way is stepping through meshio.med_io.read with pdb.runcall. Doing that, I found f['ENS_MAA']['Mesh_1']['-0000000000000000001-0000000000000000001']['MAI'][t]['FAM'][()] contains:

t n values
'SE2' 72 {0}
'TR3' 438 {-8, -7, -6}
'TE4' 839 {-9}
gdmcbain commented 5 years ago

The corresponding names for the three families of triangles and one family of tetrahedra seem to be encoded in

f['FAS']['Mesh_1']['ELEME'].keys()

<KeysViewHDF5 ['FAM_-6top', 'FAM-7bottom', 'FAM-8walls', 'FAM-9_tungsten']>

and also in, e.g.

f['FAS']['Mesh_1']['ELEME']['FAM_-6_top']['GRO']['NOM'][()]

array([[116, 111, 112, 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, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=int8)

since

 ''.join(map(chr, [116, 111, 112]))

is `'top'. (!)

The 72 2-node line-segments don't seem to belong to any such family.

gdmcbain commented 5 years ago
f['FAS']['Mesh_1']['ELEME']['FAM_-6_top'].attrs['NUM']

-6

RemDelaporteMathurin commented 5 years ago

@gdmcbain

I notice that Mesh_1_830.med contains 72 2-node line elements which are not in Mesh_1_830.msh. Are they significant? Both contain 270 3-coordinate points (matching to within ±1/10¹³), 438 3-node triangles (identical connectivity), and 839 4-node tetrahedra (same connectivity but with middle two points permuted).

Since the mesh can be read into FEniCS I would say they aren't significant.

The main difference is of course that the meshio.Mesh read from the MED contains no cell_data or field_data whereas that from your MSH has field_data:

{'top': array([1, 2]), 'bottom': array([2, 2]), 'walls': array([3, 2]), 'tungsten': array([4, 3])}

and cell_data[t]['gmsh:physical'] for t = 'triangle' (438 values in {1, 2, 3}) and tetra (839 values, all 4).

I don't really know the difference between cell_data and field_data, but the data I'm interested in are the geometrical markers for top, bottom, walls, tungsten, respectively : 1,2,3,4.

f['FAS']['Mesh_1']['ELEME']['FAM_-6_top'].attrs['NUM']

In MED format, several groups (top, bottom, walls, etc.) can belong to a family (eg. FAM_-6_top). In this specific case, you should be able to find a key named GRO containing the attribute NOM.

nschloe commented 5 years ago

I don't really know the difference between cell_data and field_data

field_data is VTK nomenclature. It's data that not associated with other mesh entities, like creation date, description strings, meta things like that. The name field_data is perhaps not well chosen.

RemDelaporteMathurin commented 5 years ago

entities should be stored in cell_data then?

nschloe commented 5 years ago

Depends on what you mean by "entities". Cell data is data associated with cells.

RemDelaporteMathurin commented 5 years ago

Sorry for that, by entity I mean this information

{'top': array([1, 2]), 'bottom': array([2, 2]), 'walls': array([3, 2]), 'tungsten': array([4, 3])}

assigned for every cell and facet. In other words, "mark" the cells and facets with their topological information. Above, the all the facets belonging to surface top are marked 1 and all the cells belonging to volume tungstenare marked 4. Is this more clear ?

nschloe commented 5 years ago

Is this more clear ?

Not to me.

assigned for every cell and facet. I

Every cell and facet has all the information? No. What does 'top': array([1, 2]) mean? Cells 1 and 2 are on the top? Or facets 1 and 2?

RemDelaporteMathurin commented 5 years ago

Every cell and facet has all the information? No. What does 'top': array([1, 2]) mean? Cells 1 and 2 are on the top? Or facets 1 and 2?

Should have said "assigned for each cell", sorry. 'top': array([1, 2]) means that the group "top" is of dimension 2, and marked with the int 1.

nschloe commented 5 years ago

'top': array([1, 2]) means that the group "top" is of dimension 2, and marked with the int 1.

Ah, okay. So that's not really cell data. I think we should find a better name for the field_data, perhaps meta_data or extra_data or something. @gdmcbain What do you think?

tianyikillua commented 5 years ago

BTW I've never understood field_data 😄

nschloe commented 5 years ago

Yeah, sorry guys! Here's the extract from the VTK manual:

vtk

tianyikillua commented 5 years ago

Okay...as you mentioned it I've indeed seen this thing under ParaView...

Returning to reading/writing point/cells sets from a MED file. The first is easy, we just need to create a new dict point_set similar to point_data. Noting that in general a same node can indeed belong to differents node sets, so point set name -> point index is appropriate.

mesh.point_set = {"set A": [0, 5, 19], 3: [2, 5, 9]}

so the point set "set A" contains indices 0, 5 and 19, while the point set 3 contains indices 2, 5 and 9, the node 5 belonging to both "set A" and 3.

tianyikillua commented 5 years ago

Concerning cell sets, I'm in favor of storing directly cell indices instead of indicator functions of each regions...as @nschloe put it in https://github.com/nschloe/meshio/issues/175#issuecomment-361647568, since "it permits cells to be a member of an arbitrary number of regions (also none at all)". With MED, I think (to be verified) this is indeed the case, where a same cell can belong to several cell sets.

Then two possibilities (https://github.com/nschloe/meshio/issues/175#issuecomment-362985744):

  1. Like cell_data: cell_type -> cell set name -> cell index
mesh.cell_set = {
    "tetra": {"set A": [0, 1, 3], 1: [1, 2]},
    "triangle": {"set B": [3, 5], 3: [0]},
    "quad": {"set C": [1, 4], 3: [0, 5, 9]},
}

where the cell set 3 is a mixed set: it contains the first triangle, and quad indices 0, 5 and 9...This set can for instance used for prescribing pressure for 3d simulations.

Also tetra index 1 belongs to two sets: "set A" and 1.

  1. Or doing cell set name -> cell_type -> cell index
mesh.cell_set = {
    "set A": {"tetra": [0, 1, 3]},
    1: {"tetra": [1, 2]},
    "set B": {"triangle": [3, 5]},
    "set C": {"quad": [1, 4]},
    3: {"triangle": [0], "quad": [0, 5, 9]}
}

Maybe the 2nd way is more readable, since in parallel with point_set we will have

PS: I omitted some hexa cells...

gdmcbain commented 5 years ago

you should be able to find a key named GRO containing the attribute NOM.

Yes. See above.

I guess that these stand for groupe (=group) and nom (=name).

nschloe commented 5 years ago

I guess that these stand for groupe (=group) and nom (=name).

There's something very particular about scientific software culture in France, and I can't say I like it.

gdmcbain commented 5 years ago

I think we should find a better name for the field_data, perhaps meta_data or extra_data or something. @gdmcbain What do you think?

I'd grown accustomed to it myself, ignoring the etymology and assuming that it was just a catch-all for any data that couldn't be associated with points or cells. (That curious assertion to the contrary in the VTK documentation notwithstanding.)

Changing the name of the attribute would break things downstream; e.g. https://github.com/kinnala/scikit-fem/blob/361a6e6475115eab8869448d2826a65a462e2e9f/skfem/mesh/mesh.py#L293-L295

though it's probably still early enough not to worry about that too much; if it grates, I think it's going to grate a long time.

I'd eschew metadata (which, were it adopted, I would set solid, without an underscore, as it's a valid word). Metadata is usually further removed than what's currently stored in field_data; metadata is usually things like the date the Mesh object was created, the name of the user, name of the host, uname, meshio version, maybe an indication of which constructor was used and what it was passed, e.g. here that meshio.read or meshio.med_io.read had been applied to a file called Mesh_1_830.med. I think, in short, that metadata is stuff that might be interesting but isn't strictly necessary. Of course, there's nothing to stop anyone putting whatever they like in field_data, including metadata like that.

On the other hand field_data contains some essential data like

{'top': array([1, 2]), 'bottom': array([2, 2]), 'walls': array([3, 2]), 'tungsten': array([4, 3])}

that really has to be passed from the upstream mesh-generator (Salome) through meshio to the finite element program so that the user can apply different boundary conditions to 'top', 'bottom' and 'walls' and appropriate coefficients to 'tungsten'.

If field_data were to be renamed, maybe just data would do. Otherwise extra_data is O. K. (other_data, misc, aux_data)?

Do any of the other formats have better names for this kind of thing?

I don't know. I don't think any of this is a compelling argument for change; but again, it is early enough to do so if desired.

gdmcbain commented 5 years ago

There's something very particular about scientific software culture in France, and I can't say I like it.

Ha, well. MEDIT .mesh on the other hand is quite nice—lean and elegant. I was going to say that on looking into MED for the first time yesterday, it looked a lot like a mesh format designed by an atomic energy agency. Very heavy and bureaucratic, like Exodus II. But I digress.

gdmcbain commented 5 years ago

Concerning cell sets, I'm in favor of storing directly cell indices instead of indicator functions of each regions...as @nschloe put it in #175 (comment), since "it permits cells to be a member of an arbitrary number of regions (also none at all)". With MED, I think (to be verified) this is indeed the case, where a same cell can belong to several cell sets.

No. MED, at least in the example provided, gives a mapping elementregion; e.g.

f['ENS_MAA']['Mesh_1']['-0000000000000000001-0000000000000000001']['MAI']['TR3']['FAM'][()]

is an array of ints, one per triangle, being one of the three labels −6, −7, −8, of the regions 'top', 'bottom' and 'walls'.

tianyikillua commented 5 years ago

Another way is to think field_data as a format-specific data: all data that can not be contained in a point_data or a cell_data, specific to a particular format.

For MED, we could store in field_data the Gauss point coordinates for instance, with which we would be able to read and write MED files with ELGA fields (data defined at Gauss points). Currently meshio can only read ELGA.

I'm okay with its name...

gdmcbain commented 5 years ago

As in #175 and #290, I do like the idea of mappings regionelement, and think that they should be supported wherever they're found, but it seems that that's not very many places. There seems to be partial support for it in Exodus II. Actually Gmsh's MSH4 supports it, but when I adapted meshio's machinery for reading MSH2 to MSH4, I coerced the data backwards to make it compatible with downstream products using pygmsh. :( Might have missed an opportunity there. Also the current meshio support for MSH4 might fail for meshes making more use of the new entity concept.

What I think might be the way for meshio to go is:

I think that that might have been aired in #175 or #290 or somewhere, but it's been a while.

gdmcbain commented 5 years ago

Another way is to think field_data as a format-specific data: all data that can not be contained in a point_data or a cell_data, specific to a particular format.

No. field_data is (at least currently) the only place to store mappings regionelement. Where else would one put them? They don't have the right shape for cell_data. This mapping should not be thought of as format-specific if meshio is to have a chance of converting between formats, as required in this issue.

Of course any truly format-specific data can also be put into field_data, as it is a catch-all. Do you mean to provide a new place to store mappings regionelement?

Perhaps somewhat unfortunately, meshio has ended up with (#175, #199) the very format-specific sounding name "gmsh:physical" as its way of storing elementregion mappings. This is currently used also in

https://github.com/nschloe/meshio/blob/master/meshio/mdpa_io.py#L195

and probably should have been in

https://github.com/nschloe/meshio/blob/master/meshio/medit_io.py#L107

but where I baulked at calling something in MEDIT "gmsh" and introduced another format-specific sounding name for something that's really synonymous… That's going to be annoying if anyone wants to attempt to convert to or from MEDIT .mesh or Gmsh MSH; it means checking for either "gmsh:physical" or "medit:ref".

And here in #347 we have in MED a third example of an elementregion mapping (currently ignored by meshio.read). Should this also be called "gmsh:physical"? This does have to be faced if one is convert MED to MSH, preserving tagged subsets of elements #290 like 'top', 'bottom', 'walls', and 'tungsten'.

Perhaps this is the opportunity to replace the keys "gmsh:physical" and "medit:ref" in cell_data[cell_type] with a single more format-neutral synonym before implementing reading "FAM" from MED.

gdmcbain commented 5 years ago

Returning to reading/writing point/cells sets from a MED file. The first is easy, we just need to create a new dict point_set similar to point_data. Noting that in general a same node can indeed belong to differents node sets, so point set name -> point index is appropriate.

mesh.point_set = {"set A": [0, 5, 19], 3: [2, 5, 9]}

so the point set "set A" contains indices 0, 5 and 19, while the point set 3 contains indices 2, 5 and 9, the node 5 belonging to both "set A" and 3.

@tianyikillua , is this kind of thing present in the examples of MED attached above? I can't find it, whereas I do see "FAM" and "NOM".

gdmcbain commented 5 years ago

I notice that Mesh_1_830.med contains 72 2-node line elements which are not in Mesh_1_830.msh. Are they significant?

Since the mesh can be read into FEniCS I would say they aren't significant.

This is tricky. So Gmsh appears to decide to ignore these "SE2" lines. Should meshio.med_io.read do the same thing? I'm thinking in the context of testing conversion #348, so that

should produce equal or equivalent meshio.Mesh objects.

One way to do this would be to discard the cells of mai[med_type] unless mai[med_type]["FAM"][()].any() (it's 72 zeros for "SE2"); however, this would break the existing write_read tests because meshio.med_io.write doesn't write les familles at all. Is the current writer correct? Are families optional? Is it only if there are some families specified that elements of family zero should be discarded?

RemDelaporteMathurin commented 5 years ago

Is the current writer correct? Are families optional? Is it only if there are some families specified that elements of family zero should be discarded?

As I showed above, even though the meshio med writer ignores the families, the mesh is readable all the same (except for les familles et les groupes). If that's your question.

RemDelaporteMathurin commented 5 years ago

if region → element is required, reverse the previous, i.e. try field_data, then `cell_data, then skip.

@gdmcbain I think this would be more flexible as we'd "only" have to create a new set whereas a mapping element → region could be a bit trickier as elements can belong to several GROUPS (not Families).

is this kind of thing present in the examples of MED attached above? I can't find it, whereas I do see "FAM" and "NOM".

This shouldn't be present in the MED files I produced cause I had no elements belonging to different groups. I should be able to produce some new examples though.

tianyikillua commented 5 years ago

Using salome I created a more complicated (way more...) mesh of a cylinder, with lots of node sets and cell sets.

The cylinder is composed of

Salome file

Cylinder_salome.zip

MED file

Cylinder.zip

It can be successfully read by meshio, without all node/cell sets.

Node sets

Edge set

Face sets

Volume sets

tianyikillua commented 5 years ago

INDEED MED IS VERY STRANGE!!! @RemiTheWarrior is from this famous atomic energy agency 😄

  1. Same nodes/elements belonging to different sets are indeed possible (see the previous MED file)
  2. These sets are organized in the following way in the HDF5 file, where internally some set intersections are created for the node/elements mentioned in 1
  3. Then each node/element receives a unique set id

tianyikillua commented 5 years ago

Why allowing same elements/nodes belonging to different sets? Maybe Sartre can answer this.

Anyway due to superposition, the engineer may prefer to apply

  1. The atmosphere pressure on Top and Down
  2. And some internal pressure on Top

Not a good example, but you see what I mean...

RemDelaporteMathurin commented 5 years ago

@RemiTheWarrior is from this famous atomic energy agency

@tianyikillua I am 😊 But I kinda agree that there's a heavy bureaucracy culture... 😉

Not a good example, but you see what I mean...

I definitely see what you mean. I should try to do the a priori sucessful conversion I did with overlapping groups to see the outputs. Maybe GMSH can't handle that

tianyikillua commented 5 years ago

Returning to reading/writing point/cells sets from a MED file. The first is easy, we just need to create a new dict point_set similar to point_data. Noting that in general a same node can indeed belong to differents node sets, so point set name -> point index is appropriate.

mesh.point_set = {"set A": [0, 5, 19], 3: [2, 5, 9]}

so the point set "set A" contains indices 0, 5 and 19, while the point set 3 contains indices 2, 5 and 9, the node 5 belonging to both "set A" and 3.

@tianyikillua , is this kind of thing present in the examples of MED attached above? I can't find it, whereas I do see "FAM" and "NOM".

I created a more complicated mesh where you can find node sets, under

FAS\[mesh name]\NOEUD

noeud simply means nodes...

tianyikillua commented 5 years ago

Is the current writer correct? Are families optional? Is it only if there are some families specified that elements of family zero should be discarded?

As I showed above, even though the meshio med writer ignores the families, the mesh is readable all the same (except for les familles et les groupes). If that's your question.

Yes the family FAM (node/element -> region id) for nodes NOE and elements under MAI is optional.

Also, the FAS indicating these regions is also optional, except the FAMILLE_ZERO.

Using meshio I read and rewrite the same cylinder file, resulting in a MED without any groups under salome.

Cylinder_meshio.zip

We can regard salome as a reference implementation of MED. The MED without groups can be read by gmsh, while the original MED file with lots of groups (and intersections of groups) can NOT be read by gmsh.

RemDelaporteMathurin commented 5 years ago

the original MED file with lots of groups (and unions of groups) can NOT be read by gmsh.

If you're using the very last version of SALOME (9.2.2) there's a known bug which automatically exports .med in version 4.0. MED 4.0 cannot be read by GMSH whether there are groups or not. So far.

I used an older version of SALOME in order to produce the example above. I'm not sure that the number of groups is an issue as long as they don't overlapp one another.

Moreover, in the conversion above, I noticed that if all the cells aren't in a group, GMSH just ignores this cell, resulting in a hollow mesh.

tianyikillua commented 5 years ago

So for MED, there are two possibilities of implementing cell sets, node sets being trivial.

  1. Using user-defined sets and internally created intersections sets for overlapped sets to define a element -> region id (currently done in meshio via gmsh_physical as a cell_data and in MED file internally, through FAM)
  2. Go deeper into this element -> region id by decomposing possible intersections sets, and associate several groups to the same node/element if it is the case. Construct a mesh.cell_set.
tianyikillua commented 5 years ago

the original MED file with lots of groups (and intersections of groups) can NOT be read by gmsh.

If you're using the very last version of SALOME (9.2.2) there's a known bug which automatically exports .med in version 4.0. MED 4.0 cannot be read by GMSH whether there are groups or not. So far.

I used an older version of SALOME in order to produce the example above. I'm not sure that the number of groups is an issue as long as they don't overlapp one another.

Moreover, in the conversion above, I noticed that if all the cells aren't in a group, GMSH just ignores this cell, resulting in a hollow mesh.

Oh thanks! You were right. Below is the same cylinder MED file, where the version is changed from 4.0.0 to 3.0.0, by modifying manually INFOS_GENERALES.

Cylinder_3_0_0.zip

This file can be loaded by gmsh (it and knows treating overlapping sets and elements belonging to different groups, I think).

Using gmsh to convert MED to msh, and read this msh using meshio gives the above figure.

{'Top circle': array([4, 1]),
 'Top': array([5, 2]),
 'Top and down': array([6, 2]),
 'B': array([1, 3]),
 'A': array([2, 3]),
 'C': array([3, 3])}

which is not useful since gmsh:physical is all 0. Where is the field_data for gmsh:geometrical?

RemDelaporteMathurin commented 5 years ago
  • gmsh:physical is all 0, since we don't know how to treat elements belonging to several groups (good job)

@tianyikillua So we're able to read those belonging to only one group then ?

tianyikillua commented 5 years ago
  • gmsh:physical is all 0, since we don't know how to treat elements belonging to several groups (good job)

@tianyikillua So we're able to read those belonging to only one group then ?

Cylinder_no_overlap_groups_3_0_0.zip

Sadly no...still 0 everywhere. In the msh file the elements refer to the gmsh:geometrical region, which doesn't match gmsh:physical contained in field_data.