Manipuler les DataArray
-----------------------
Les DataArrays (``DataArrayInt`` et ``DataArrayDouble``) sont utilisĂŠs dans MEDCoupling pour stocker des valeurs sous
forme de tableaux contigus en mĂŠmoire. Les valeurs sont groupĂŠes par tuples, et chaque tuple a le mĂŞme nombre de composantes.
Ils sont Ă la base de beaucoup de traitements rĂŠalisĂŠs dans MEDCoupling. Il est ainsi important de bien savoir les manipuler.
Les ``DataArrayDouble`` sont souvent utilisĂŠs pour la manipulation directe des valeurs d'un champ comme on le verra plus tard.
Les ``DataArrayInt`` eux sont utilisĂŠs pour toutes les fonctionnalitĂŠs travaillant avec des identifiants de
cellules et/ou de points.
Le but de l'exercice
~~~~~~~~~~~~~~~~~~~~
Le but ici est de crĂŠer les coordonnĂŠes de 7 hexagones rĂŠguliers (tous inscrits dans des cercles de rayon 3m) en dimension 2.
La première composante du tableau de coordonnÊes s'appelera X avec l'unitÊ "m" (mètre) et la 2ème composante s'appelera "Y" avec la même unitÊ.
On pourrait directement calculer les coordonnĂŠes de l'ensemble des points requis avec un peu de trigonomĂŠtrie, mais afin de travailler un peu avec l'API, on fait le choix de construire les 7 hexagones Ă partir d'un seul hexagone rĂŠgulier centrĂŠ en [3.4; 4.4].
Autour de cet hexagone rÊgulier central, on crÊe 6 copies translatÊes et chaque copie partagera exactement un bord (edge) avec le motif initial. Ensuite on fusionne les noeuds (tuples) communs. Ceci nous permettra de manipuler les indirections et les mÊthodes d'indexing très usitÊes dans les maillages non structurÊs.
.. image:: images/DataArrayDouble_1.jpg
:scale: 50
Les points traitĂŠs ici :
* CrĂŠer une instance de ``DataArrayDouble``
* Afficher une instance de ``DataArrayDouble`` et invoquer la mĂŠthode ``getValue()`` pour la convertir en liste
* Utiliser les notations pratiques ``da[:,:]`` ...
* Apprendre la renumĂŠrotation (convention "old-2-new")
* Invoquer des services tels que ``findCommonTuples()``
DĂŠbut de l'implĂŠmentation
~~~~~~~~~~~~~~~~~~~~~~~~~
Pour commencer l'exercice importer le module Python ``MEDCoupling`` et l'aliaser avec ``mc`` (ca nous ĂŠvitera des noms trop longs). Importer aussi le module ``math``. ::
import MEDCoupling as mc
import math
On rappelle que toutes les mĂŠthodes statiques du module commencent par une majuscule.
Avec ces imports sont disponibles :
* toutes les classes de MEDCoupling
* tous les ĂŠnumĂŠrations (par exemple, les types de cellules standard: ``mc.ON_CELLS``, ``mc.ON_NODES``, ``mc.ONE_TIME``...)
* toutes les mĂŠthodes statiques
CrĂŠer une instance de DataArrayDouble contenant 6 tuples
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Le but ici est de crĂŠer un ``DataArrayDouble`` contenant les coordonnĂŠes d'un seul hexagone rĂŠgulier. ::
d = mc.DataArrayDouble(6,2)
Ceci est ĂŠquivalent Ă ::
d = mc.DataArrayDouble()
d.alloc(6,2)
Ceci est aussi ĂŠquivalent Ă ::
d = mc.DataArrayDouble(12)
d.rearrange(2)
Notons enfin que l'on peut aussi directement construire un ``DataArray`` Ă partir d'une liste Python. Par dĂŠfaut le tableau
n'a qu'une seule composante. ::
d_example = mc.DataArrayDouble([0.0,1.0,2.5])
print d_example
.. note:: Le tableau ``d`` contient maintenant 12 valeurs groupĂŠes en 6 tuples contenant chacun 2 composantes.
Les valeurs dans ``d`` ne sont pas encore assignĂŠes.
Initialiser une instance de DataArrayDouble
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Assigner la valeur 3.0 (le rayon) à la première composante de chacun des tuples. La syntaxe ressemble fortement
à celle de NumPy. On peut par exemple assigner d'un coup les tuples 1 à 3 (inclus), sur la première composante avec la valeur
3.0 ::
d[1:4,0] = 3.
ou directement l'intÊgralitÊ de la première composante ::
d[:,0] = 3.
Initialiser la 2ème composante de chaque tuple i avec la valeur i. ::
d[:,1] = range(6)
Multiplier la seconde composante de chacun des tuples par pi/3. ::
d[:,1] *= math.pi/3.
.. note:: ``d`` contient dĂŠsormais les coordonnĂŠes polaires des noeuds de notre hexagone rĂŠgulier centrĂŠ en 0,0 pour le moment.
Convertir ``d`` de polaire Ă cartĂŠsien en invoquant la mĂŠthode ``fromPolarToCart()`` et re-mettre le rĂŠsultat dans ``d``. ::
d = d.fromPolarToCart()
.. note:: ``fromPolarToCart()`` gÊnère une nouvelle instance, nous avons donc perdu le ``d`` initial
Assigner les informations textuelles correctes sur les 2 composantes de ``d`` : ::
d.setInfoOnComponents(["X [m]","Y [m]"])
.. note:: Cela n'est pas indispensable pour cet exercise, mais d'autres fonctions plus avancĂŠes nĂŠcessitent cette information.
Afficher ``d`` tel quel. ::
print d
Afficher juste les valeurs sous forme d'une liste python. ::
print d.getValues()
VÊrifier que pour chaque tuple dÊsormais dans ``d``, sa norme (mÊthode ``magnitude()``) est bien Êgale à 3.0, à 1.e-12 près (mÊthode ``isUniform()``) ::
print "Uniform array?", d.magnitude().isUniform(3.,1e-12)
Duplication et agrĂŠgation
~~~~~~~~~~~~~~~~~~~~~~~~~
On construit maintenant la liste ``translationToPerform``, qui contient une liste de vecteurs chacun de taille 2. Cette liste de taille 7 (7 hexagones) contient la translation Ă opĂŠrer pour produire chacun des hexagones.
Faites nous confiance sur la trigonomĂŠtrie, vous pouvez copier directement les deux lignes suivantes ::
radius = 3.
translationToPerform = [[0.,0.],[3./2.*radius,-radius*math.sqrt(3.)/2],[3./2.*radius,radius*math.sqrt(3.)/2],[0.,radius*math.sqrt(3.)],[-3./2.*radius,radius*math.sqrt(3.)/2],[-3./2.*radius,-radius*math.sqrt(3.)/2],[0.,-radius*math.sqrt(3.)]]
CrĂŠer les 7 copies de ``d`` et opĂŠrer la "translation" correspondante. ::
ds = len(translationToPerform)*[None]
for pos,t in enumerate(translationToPerform):
ds[pos] = d[:] # Perform a deep copy of d and place it at position 'pos' in ds
ds[pos] += t # Adding a vector to a set of coordinates does a translation. t could have been a DataArrayDouble too.
pass
.. note:: Le ``pass`` Ă la fin de la boucle ``for`` n'est pas indispensable mais aide certains ĂŠditeurs Ă indenter le code.
Une autre façon de faire un peu plus compacte (pour les amoureux des *one-liner*) : ::
ds = [d + translationToPerform[i] for i in xrange(len(translationToPerform))]
AgrĂŠgation de tableaux
~~~~~~~~~~~~~~~~~~~~~~
A partir de la liste d'instances de DataArrayDouble ``ds`` construire le DataArrayDouble ``d2`` rĂŠsultat de l'*agrĂŠgation* des instances les unes Ă la suite des autres. ::
d2 = mc.DataArrayDouble.Aggregate(ds)
``d2`` contient dĂŠsormais l'ensemble des tuples (6*7 de 2 composantes chacun) des
instances contenues dans ``ds``, en respectant l'ordre dans ``ds``. Cela parait ĂŠvident, mais
l'agrĂŠgation de maillages et de champs respecte exactement le mĂŞme principe pour
faciliter l'accès et le repÊrage des donnÊes. C'est par exemple une diffÊrence essentielle avec le
modèle MED fichier comme on le verra plus tard.
.. note:: La mĂŠthode permettant d'agrĂŠger par composante (c'est-Ă -dire de concatĂŠner des tableaux
colonne par colonne, plutĂ´t que par tuples) s'appelle ``Meld()``.
Trouver les tuples ĂŠgaux
~~~~~~~~~~~~~~~~~~~~~~~~
La variable ``d2`` contient 42 tuples mais certains tuples apparaissent plusieurs fois.
Pour trouver les tuples Êgaux à 1e-12 près (prÊcision absolue) invoquer ``findCommonTuples()``.
Utiliser ``help(mc.DataArrayDouble.findCommonTuples)`` pour en connaitre l'interface. Stocker le retour de la fonction dans
``c`` et ``cI`` ::
oldNbOfTuples = d2.getNumberOfTuples()
c,cI = d2.findCommonTuples(1e-12)
On a ainsi rĂŠcupĂŠrĂŠ dans ``c`` l'ensemble des m=12 groupes de noeuds communs accollĂŠs. ``cI`` contient les index pour repĂŠrer les identifiants de points dans ``c`` pour tout groupe
``i`` dans [0,12). Ainsi les identifiants de tuples du groupe ``i`` commencent Ă l'index ``cI[i]`` et finissent Ă l'index ``cI[i+1]``.
La mÊthode ``findCommonTuples()`` retourne ainsi 2 paramètres: un tableau contenant la liste des tuples communs
et un tableau d'index qui permet de naviguer dans le premier tableau.
Il s'agit d'une forme de retour très classique dans MEDCoupling, appelÊe *indirect indexing*. Cela apparaÎt souvent dans la manipulation des
maillages non structurĂŠs. Cette reprĂŠsentation est rappelĂŠe sur l'image ci-dessous, oĂš le premier tableau et en haut,
et le deuxième tableau permettant de la parcourir en bas:
.. image:: images/IndirectIndex.jpg
:scale: 50
.. note:: Le dernier ĂŠlĂŠment de ``cI`` pointe en dehors du tableau ``c``. Ce dernier index est toujours prĂŠsent
et permet de s'assurer que des traitements tels que les *slices* prÊsentÊs juste après, sont toujours valables,
sans avoir besoin de particulariser le dernier groupe.
.. _indirect-index-exo:
Manipuler le format "indirect index"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Le nombre de tuples communs à 1e-12 près est donc Êgal à ``len(cI)-1``, c'est-à -dire 12 dans notre cas.
RĂŠcupĂŠrer la liste des identifiants de tuples du groupe 0 et mettre le rĂŠsultat dans la variable ``tmp``.
Afficher ``tmp``. ::
tmp = c[cI[0]:cI[0+1]]
print tmp
VĂŠrifier, en l'affichant, que pour tous les identifiants de tuples dans ``tmp``, leurs tuples sont bien ĂŠgaux dans ``d2``. ::
print d2[tmp]
.. note:: On voit que le tuple (3.,0.) à 1e-12 près est rÊpÊtÊ 3 fois et ``tmp`` donne les positions respectives de
ces 3 rĂŠpĂŠtitions.
Maintenant on va dĂŠduire des variables ``oldNbOfTuples``, ``c`` et ``cI`` le nombre de tuples effectivement diffĂŠrents dans d2.
Pour ce faire, nous allons trouver le nombre de tuples doublons dans ``d2`` et soustraire le rĂŠsultat de ``oldNbOfTuples``.
Pour connaĂŽtre le nombre de doublons, invoquer ``DataArrayInt.deltaShiftIndex`` qui retourne pour chaque groupe sa taille.
Mettre le rĂŠsultat dans ``a``. ::
a = cI.deltaShiftIndex()
DĂŠduire de ``a`` le nombre de tuples doublons dans ``d2`` par groupe et mettre le rĂŠsultat dans ``b``. ::
b = a-1
Enfin on peut trouver le nouveau nombre de tuples grâce à ``b`` et à ``oldNbOfTuples``. Mettre le rÊsultat dans ``myNewNbOfTuples``. ::
myNewNbOfTuples = oldNbOfTuples - sum(b.getValues())
Construire un tableau "old-2-new"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Nous allons maintenant exploiter cette information pour extraire un seul
reprĂŠsentant dans chaque groupe de points dupliquĂŠs.
Les deux tableaux ``c`` et ``cI`` dĂŠfinissent une surjection d'un espace de dĂŠpart Ă 42 (``oldNbOfTuples``) tuples X
vers un espace Ă 24 (``myNewNbOfTuples``) tuples Y.
.. image:: images/SurjectionDataArray.png
L'autre manière de dÊfinir cette surjection (sans perte d'information) est de la reprÊsenter par un tableau "old-2-new".
Ce mode de stockage prend la forme d'un DataArrayInt ``o2n`` composĂŠ de Card(X) tuples (i.e. 42) Ă une composante.
Pour chaque tuple (ĂŠlĂŠment) d'index ``i`` de ``o2n``, la case ``o2n[i]`` contient le nouvel identifiant de tuple dans Y.
On va donc d'un ancien identifiant (old) vers un nouveau (new).
Nous allons construire ce tableau pour extraire un sous-ensemble des coordonnĂŠes de dĂŠpart, et ne garder que les
tuples uniques (non doublons) dans l'ensemble de dĂŠpart.
.. note:: Pour toutes les opĂŠrations de renumĂŠrotation en MEDCoupling (bijection),
le format "old-2-new" est systĂŠmatiquement utilisĂŠ.
La mĂŠthode statique ``DataArrayInt.ConvertIndexArrayToO2N()`` (nom un peu barbare, on vous l'accorde)
permet de passer du mode de stockage de cette surjection ``c``, ``cI`` au format ``o2n``.
On rÊcupère au passage card(Y) c'est-à -dire le ``newNbOfTuples``. ::
o2n, newNbOfTuples = mc.DataArrayInt.ConvertIndexArrayToO2N(oldNbOfTuples,c,cI)
print "Have I got the right number of tuples?"
print "myNewNbOfTuples = %d, newNbOfTuples = %d" % (myNewNbOfTuples, newNbOfTuples)
assert(myNewNbOfTuples == newNbOfTuples)
Nous pouvons maintenant constuire le tableau de points uniques ``d3``. A l'aide de ``o2n``
et ``newNbOfTuples``, invoquer ``DataArrayDouble.renumberAndReduce()`` sur ``d2``. ::
d3 = d2.renumberAndReduce(o2n, newNbOfTuples)
L'inconvĂŠnient de cette mĂŠthode c'est que finalement on ne connait pas pour chaque groupe de tuple communs dans
d2 quel identifiant a ĂŠtĂŠ utilisĂŠ.
Par exemple pour le groupe 0 on sait que les tuples 0, 8 et 16 (tmp.getValues()) sont tous ĂŠgaux, et on ne sait
pas si 0, 8 ou 16 a ĂŠtĂŠ utilisĂŠ pour remplir ``d3``.
Si l'on souhaite expliciter ce choix, on peut passer en format "new-2-old". Ce mode de stockage prend la forme d'un
``DataArrayInt`` ``n2o`` composĂŠ de Card(Y)
tuples (24) Ă 1 composante. Pour chaque tuple (ĂŠlĂŠment) d'index i de ``n2o``, la case ``n2o[i]`` contient l'index du tuple qui a ĂŠtĂŠ choisi dans X.
Pour passer d'une description "old-2-new" vers "new-2-old", la mĂŠthode est ``DataArrayInt.invertArrayO2N2N2O()``.
Effectuer ce traitement sur la variable ``o2n``. ::
n2o = o2n.invertArrayO2N2N2O(newNbOfTuples)
A l'aide de ``n2o`` on peut construire un ``d3_bis`` Ă partir de ``d2``, et qui contient la mĂŞme chose que le ``d3`` prĂŠcĂŠdent. ::
d3_bis = d2[n2o]
print "Are d3 and d3_bis equal ? %s" % (str(d3.isEqual(d3_bis, 1e-12)))
Translater tous les tuples
~~~~~~~~~~~~~~~~~~~~~~~~~~
Tous les tuples (ou nodes) sont Ă translater du vecteur [3.3,4.4] afin de recentrer toute la figure en ce point. ::
d3 += [3.3,4.4]
Constuire un maillage non strucuturĂŠ
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
On chercher maintenant Ă crĂŠer le maillage final montrĂŠ dans la figure. Nous avons dĂŠjĂ construit le tableau
de coordonnĂŠes, il nous reste les cellules Ă crĂŠer.
CrĂŠer un maillage non structurĂŠ ``m`` avec les coordonnĂŠes ``d3``. Le maillage``m`` a une mesh-dimension 2 ::
m = mc.MEDCouplingUMesh("My7hexagons",2)
m.setCoords(d3)
print "Mesh dimension is", m.getMeshDimension()
print "Spatial dimension is", m.getCoords().getNumberOfComponents()
Maintenant, allouer le nombre de cellules avec (un majorant du) nombre attendu de cellules. ::
m.allocateCells(7)
Enfin grâce à ``o2n`` on a la *connectivitÊ* (i.e. la liste des points formant un hexagone)
des 7 hexagones utilisant les coordonnĂŠes ``d3``. ::
for i in xrange(7):
cell_connec = o2n[6*i:6*(i+1)]
m.insertNextCell(mc.NORM_POLYGON, cell_connec.getValues())
pass
VĂŠrifier que ``m`` est correct et ne contient pas d'anomalie. ::
m.checkConsistencyLight()
.. note:: Il est toujours une bonne idÊe d'appeler cette mÊthode après la construction "from scratch" d'un maillage.
Cela assure qu'il n'y a pas de gros "couacs" dans la connectivitĂŠ, etc ...
Pour vĂŠrifier *visuellment* que ``m`` est correct, l'ĂŠcrire dans un fichier "My7hexagons.vtu" et le visualiser dans ParaViS. ::
m.writeVTK("My7hexagons.vtu")
.. note:: On a ĂŠcrit ici dans un fichier VTU et non MED, car MEDCoupling n'inclut pas par dĂŠfaut les services de MED fichier.
Bien que l'on ĂŠcrive au format VTK (\*.vtu), MEDCoupling ne dĂŠpend pas de VTK.
Solution
~~~~~~~~
:ref:`python_testMEDCouplingdataarray1_solution`