[VTK] Écrire un fichier VTI

Dans cette partie, on s'intéressera au maillages dit structurés et aux différentes manières de les sauvegarder dans des fichiers vtk. En fonction des caractéristiques du maillage, ces fichiers auront l'exetension .vti, .vtr ou .vts.

Après une petite introduction sur les maillages structurés, passons à la pratique.

Le format VTI est le format xml de vtk pour sauvegarder des données évaluées dans un maillage structuré cartésien appelé ImageData ou

S’il te plaît dessine moi un ImageData

Ce maillage est défini par :

  • son origine (x0, y0, z0),
  • les dimensions d’une cellule (dx, dy, dz)
  • et les dimensions du maillage (nx, ny, nz). Les dimensions correspondent aux nombres de cellules dans chaque direction.
<?xml version="1.0"?>
<VTKFile type="ImagData" version="0.1" byte_order="LittleEndian" ...>
  <ImageData WholeExtent="x1 x2 y1 y2 z1 z2"
   Origin="x0 y0 z0" Spacing="dx dy dz">
   <Piece Extent="x1 x2 y1 y2 z1 z2">
      <PointData>...</PointData>
      <CellData>...</CellData>
   </Piece>
   </ImageData>
</VTKFile>

La balise VTKFile a attribut type égale à ImageData. L’attribut byte_order a deux valeurs possible LittleEndian ou BigEndian. Cette valeur dépend de votre système et vous pouvez la récupérer comme suit :

print(sys.byteorder)

Il faut une balise ImageData dont les attributs Origin (origine du maillage), Spacing (dimension d’une cellule) et WholeExtent (dimension du maillage) suffisent à décrire le maillage. Et enfin une balise Piece dont l’attribut Extent sera la valeur de WholeExtent pour le type de vti que nous écrirons; un vti qui ne décrit pas un domaine découpé en sous-domaine pour la parallélisation. Ainsi la plupart du temps “x1 x2 y1 y2 z1 z2” vaut “0 nx 0 ny 0 nz”. Enfin les données à représenter sont entre les balises PointData et CellData. Ces balises n’apparaissent qu’une seule fois même si vous avez plusieurs données de type pointdata ou celldata.

Écrire les données

Remarquez qu’un tableau contenant des PointData doit contenir exactement (nx+1)(ny+1)(nz+1) éléments et nx*ny*nz pour un tableau contenant des CellData. Une données est écrite entre les balises DataArray et écrite en ascii ou en binary: Il y a autant de DataArray que de données entre les balises PointData et CellData.

Prenons l’exemple d’un tableau donné aux 10 noeuds de mon maillage:

<DataArray type="Int32" Name="temperature" format="ascii">
0 1 2 3 4 5 6 7 8 9
</DataArray>
<DataArray type="Float32" Name="cell_normals" NumberOfComponents="3" format="ascii">
0 0 -1 0 0 1 0 -1 0 0 1 0 -1 0 0 1 0 0 0 -1 0 0 1 0 -1 0 0 1 0 0
</DataArray>

Voici l’exemple d’un tableau écrit sous les 2 formes. Il s’agit exactement du même tableau. Le tableau

<DataArray Name="density" NumberOfComponents="1" type="Float64" format="ascii">
0 1 2 3 4 5 6 7
</DataArray>
<DataArray Name="dens" NumberOfComponents="1" type="Float64" format="binary">
QAAAAAAAAAA=AAAAAAAAAAAAAAAAAADwPwAAAAAAAABAAAAAAAAACEAAAAAAAAAQQAAAAAAAABRAAAAAAAAAGEAAAAAAAAAcQA==
</DataArray>
</CellData>

Les attribut de DataArray sont:

  • type — le type de la donnée Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64.
  • Name — le nom de la donnée.
  • NumberOfComponents — Est-ce une donnée 1D, 2D ou 3D.
  • format — comment est sauvegardée la donnée dans le fichier; “ascii”, “binary”, ou “appended”.
  • offset — si le format est “appended”, il faut spécifier la position de la première valeur du tableau à partir du symbole “_“.

Du choix de attribut format dépend de la manière et où sera écrit les données :

  • format=”ascii” — les valeurs sont écrites directement entre les balises DataArray et séparées par des espaces.
  • format=”binary” — les valeurs sont encodées en base64 et écrites sans séparation directement entre les balises DataArray. L’ordre byte-order doit respecter celui spécifier dans la balise VTKFILE.
  • format=”appended” — toute les données avec cette valeur sont écrites en binaire et à la chaîne dans une unique balise, AppendedData. L’attribut offset permet de savoir où commence la donnée qui nous intéresse.

Attention : en binaire avant d’écrire votre tableau, il faut écrire la taille (UInt64) en byte de votre tableau. Ainsi, on retrouve dans l’exemple précédent QAAAAAAAAAA=AAAAAAAAAAAAAAAA… QAAAAAAAAAA= étant la taille du tableau de float64 soit 8*8

print( struct.unpack("<Q",base64.b64decode(b"QAAAAAAAAAA=")) )

_QMwEAAAAAAAAA…

Si vous choisissez le format appended, avant d’écrire le premier tableau, il faut mettle symbole “_” underscore. Ce symbole ne fait pas parti des données mais est obligatoire. Chaque donnée des DataArray est sauvegardé en continue (l’un après les autres) sans élément séparateur.

Écrire un fichier vti en Python

import numpy as np
import struct, base64
xO, yO, zO = [1., -2., 0.]
nx, ny, nz = 2, 2, 2
dx, dy, dz = 1., 2., 1.

npoints = (nx+1)*(ny+1)*(nz+1)
ncells  = nx*ny*nz

temp    = np.arange(npoints, dtype=np.float64)
density = np.random.rand(ncells)
dens    = np.arange(ncells)

def block_size(array):
    return base64.urlsafe_b64encode(struct.pack("<Q", 8*array.size)) # UInt64

fname = "image.vti"
with open(fname, "w") as txt:
    txt.write('<?xml version="1.0"?>\n')
    txt.write(f'<VTKFile type="ImageData" version="1.0" byte_order="LittleEndian" header_type="UInt64">\n')
    #txt.write('<VTKFile type="ImageData" version="1.0" byte_order="LittleEndian">\n')
    txt.write(f' <ImageData WholeExtent="0 {nx} 0 {ny} 0 {nz}" Origin="{xO} {yO} {zO}" Spacing="{dx} {dy} {dz}">\n')
    txt.write(f'  <Piece Extent="0 {nx} 0 {ny} 0 {nz}">\n')
    txt.write('  <PointData>\n')
    txt.write('  <DataArray Name="temp2" NumberOfComponents="1" type="Float64" format="asci">\n')
    for t in temp:
        print(t, end=" ", file=txt)
    print("\n  </DataArray>", file=txt)
with open(fname, "ba") as txt:
    txt.write(b'  <DataArray Name="temp" NumberOfComponents="1" type="Float64" format="binary">\n')
    txt.write(block_size(temp))
    txt.write(base64.b64encode(struct.pack(f"<{temp.size}d", *temp)))
    txt.write(b"\n  </DataArray>\n")
with open(fname, "a") as txt:
    txt.write('  </PointData>\n')
    txt.write('  <CellData>\n')
    txt.write('  <DataArray Name="density" NumberOfComponents="1" type="Float64" format="ascii">\n')
    for d in dens:
        print(d, end=" ", file=txt)
    print("\n  </DataArray>", file=txt)
with open(fname, "ba") as txt:
    txt.write(b'  <DataArray Name="dens" NumberOfComponents="1" type="Float64" format="binary">\n')
    txt.write(block_size(dens))
    txt.write(base64.b64encode(struct.pack(f"<{dens.size}d", *dens)))
    txt.write(b"\n  </DataArray>\n")
    txt.write(b'  <DataArray Name="dens2" NumberOfComponents="1" type="Float64" format="binary" encoding="base64">\n')
    txt.write(block_size(dens))
    txt.write(base64.b64encode(struct.pack(f"<{dens.size}d", *dens)))
    txt.write(b"\n  </DataArray>\n")
with open(fname, "a") as txt:
    txt.write('  </CellData>\n')
    txt.write('</Piece>\n')
    txt.write('</ImageData>\n')
    txt.write('</VTKFile>\n')
<?xml version="1.0"?>
<VTKFile type="ImageData" version="1.0" byte_order="LittleEndian" header_type="UInt64">
 <ImageData WholeExtent="0 2 0 2 0 2" Origin="1.0 -2.0 0.0" Spacing="1.0 2.0 1.0">
  <Piece Extent="0 2 0 2 0 2">
  <PointData>
  <DataArray Name="temp2" NumberOfComponents="1" type="Float64" format="asci">
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0 24.0 25.0 26.0
  </DataArray>
  <DataArray Name="temp" NumberOfComponents="1" type="Float64" format="binary">
2AAAAAAAAAA=AAAAAAAAAAAAAAAAAADwPwAAAAAAAABAAAAAAAAACEAAAAAAAAAQQAAAAAAAABRAAAAAAAAAGEAAAAAAAAAcQAAAAAAAACBAAAAAAAAAIkAAAAAAAAAkQAAAAAAAACZAAAAAAAAAKEAAAAAAAAAqQAAAAAAAACxAAAAAAAAALkAAAAAAAAAwQAAAAAAAADFAAAAAAAAAMkAAAAAAAAAzQAAAAAAAADRAAAAAAAAANUAAAAAAAAA2QAAAAAAAADdAAAAAAAAAOEAAAAAAAAA5QAAAAAAAADpA
  </DataArray>
  </PointData>
  <CellData>
  <DataArray Name="density" NumberOfComponents="1" type="Float64" format="ascii">
0 1 2 3 4 5 6 7
  </DataArray>
  <DataArray Name="dens" NumberOfComponents="1" type="Float64" format="binary">
QAAAAAAAAAA=AAAAAAAAAAAAAAAAAADwPwAAAAAAAABAAAAAAAAACEAAAAAAAAAQQAAAAAAAABRAAAAAAAAAGEAAAAAAAAAcQA==
  </DataArray>
  <DataArray Name="dens2" NumberOfComponents="1" type="Float64" format="binary" encoding="base64">
QAAAAAAAAAA=AAAAAAAAAAAAAAAAAADwPwAAAAAAAABAAAAAAAAACEAAAAAAAAAQQAAAAAAAABRAAAAAAAAAGEAAAAAAAAAcQA==
  </DataArray>
  </CellData>
</Piece>
</ImageData>
</VTKFile>

Liens utiles