]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/medcoupling_dataarray1_fr.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / doc / tutorial / medcoupling_dataarray1_fr.rst
1
2 Manipuler les DataArray
3 -----------------------
4
5 Les DataArrays (``DataArrayInt`` et ``DataArrayDouble``) sont utilisés dans MEDCoupling pour stocker des valeurs sous 
6 forme de tableaux contigus en mémoire. Les valeurs sont groupées par tuples, et chaque tuple a le même nombre de composantes.
7 Ils sont à la base de beaucoup de traitements réalisés dans MEDCoupling. Il est ainsi important de bien savoir les manipuler.
8
9 Les ``DataArrayDouble`` sont souvent utilisés pour la manipulation directe des valeurs d'un champ comme on le verra plus tard.
10 Les ``DataArrayInt`` eux sont utilisés pour toutes les fonctionnalités travaillant avec des identifiants de 
11 cellules et/ou de points.
12
13 Le but de l'exercice
14 ~~~~~~~~~~~~~~~~~~~~
15
16 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.
17 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é.
18
19 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].
20 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.
21
22 .. image:: images/DataArrayDouble_1.jpg
23         :scale: 50
24
25 Les points traités ici :
26
27 * Créer une instance de ``DataArrayDouble``
28 * Afficher une instance de ``DataArrayDouble`` et invoquer la méthode ``getValue()`` pour la convertir en liste
29 * Utiliser les notations pratiques ``da[:,:]`` ...
30 * Apprendre la renumérotation (convention "old-2-new")
31 * Invoquer des services tels que ``findCommonTuples()``
32
33 Début de l'implémentation
34 ~~~~~~~~~~~~~~~~~~~~~~~~~
35
36 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``. ::
37
38         import MEDCoupling as mc
39         import math
40
41 On rappelle que toutes les méthodes statiques du module commencent par une majuscule. 
42 Avec ces imports sont disponibles :
43
44 * toutes les classes de MEDCoupling
45 * tous les énumérations (par exemple, les types de cellules standard: ``mc.ON_CELLS``, ``mc.ON_NODES``, ``mc.ONE_TIME``...)
46 * toutes les méthodes statiques
47
48 Créer une instance de DataArrayDouble contenant 6 tuples
49 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
50
51 Le but ici est de créer un ``DataArrayDouble`` contenant les coordonnées d'un seul hexagone régulier. ::
52
53         d = mc.DataArrayDouble(6,2)
54
55 Ceci est équivalent à ::
56
57         d = mc.DataArrayDouble()
58         d.alloc(6,2)
59
60 Ceci est aussi équivalent à ::
61
62         d = mc.DataArrayDouble(12)
63         d.rearrange(2)
64
65 Notons enfin que l'on peut aussi directement construire un ``DataArray`` à partir d'une liste Python. Par défaut le tableau 
66 n'a qu'une seule composante. ::
67
68         d_example = mc.DataArrayDouble([0.0,1.0,2.5])
69         print d_example 
70
71 .. note:: Le tableau ``d`` contient maintenant 12 valeurs groupées en 6 tuples contenant chacun 2 composantes.
72           Les valeurs dans ``d`` ne sont pas encore assignées.
73
74 Initialiser une instance de DataArrayDouble
75 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
76
77 Assigner la valeur 3.0 (le rayon) à la première composante de chacun des tuples. La syntaxe ressemble fortement
78 à celle de NumPy. On peut par exemple assigner d'un coup les tuples 1 à 3 (inclus), sur la première composante avec la valeur
79 3.0 ::
80
81         d[1:4,0] = 3.
82
83 ou directement l'intégralité de la première composante ::
84
85         d[:,0] = 3.
86
87 Initialiser la 2ème composante de chaque tuple i avec la valeur i. ::
88
89         d[:,1] = range(6)
90
91 Multiplier la seconde composante de chacun des tuples par pi/3. ::
92
93         d[:,1] *= math.pi/3.
94
95 .. note:: ``d`` contient désormais les coordonnées polaires des noeuds de notre hexagone régulier centré en 0,0 pour le moment.
96
97 Convertir ``d`` de polaire à cartésien en invoquant la méthode ``fromPolarToCart()`` et re-mettre le résultat dans ``d``. ::
98
99         d = d.fromPolarToCart()
100
101 .. note:: ``fromPolarToCart()`` génère une nouvelle instance, nous avons donc perdu le ``d`` initial
102
103 Assigner les informations textuelles correctes sur les 2 composantes de ``d`` : ::
104
105         d.setInfoOnComponents(["X [m]","Y [m]"])
106
107 .. note:: Cela n'est pas indispensable pour cet exercise, mais d'autres fonctions plus avancées nécessitent cette information.
108
109 Afficher ``d`` tel quel. ::
110
111         print d
112
113
114 Afficher juste les valeurs sous forme d'une liste python. ::
115
116         print d.getValues()
117
118
119 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()``) ::
120
121         print "Uniform array?", d.magnitude().isUniform(3.,1e-12)
122
123
124 Duplication et agrégation
125 ~~~~~~~~~~~~~~~~~~~~~~~~~
126
127 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.
128
129 Faites nous confiance sur la trigonométrie, vous pouvez copier directement les deux lignes suivantes ::
130
131         radius = 3.
132         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.)]]
133
134
135 Créer les 7 copies de ``d`` et opérer la "translation" correspondante.  ::
136
137         ds = len(translationToPerform)*[None]
138         for pos,t in enumerate(translationToPerform):
139                 ds[pos] = d[:]      # Perform a deep copy of d and place it at position 'pos' in ds
140                 ds[pos] += t        # Adding a vector to a set of coordinates does a translation. t could have been a DataArrayDouble too.
141                 pass
142           
143 .. note:: Le ``pass`` à la fin de la boucle ``for`` n'est pas indispensable mais aide certains éditeurs à indenter le code.
144
145 Une autre façon de faire un peu plus compacte (pour les amoureux des *one-liner*) : ::
146
147         ds = [d + translationToPerform[i] for i in xrange(len(translationToPerform))]
148
149 Agrégation de tableaux
150 ~~~~~~~~~~~~~~~~~~~~~~
151
152 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. ::
153
154         d2 = mc.DataArrayDouble.Aggregate(ds)
155
156 ``d2`` contient désormais l'ensemble des tuples (6*7 de 2 composantes chacun) des
157 instances contenues dans ``ds``, en respectant l'ordre dans ``ds``. Cela parait évident, mais
158 l'agrégation de maillages et de champs respecte exactement le même principe pour
159 faciliter l'accès et le repérage des données. C'est par exemple une différence essentielle avec le
160 modèle MED fichier comme on le verra plus tard.
161
162 .. note:: La méthode permettant d'agréger par composante (c'est-à-dire de concaténer des tableaux
163     colonne par colonne, plutôt que par tuples) s'appelle ``Meld()``.
164
165 Trouver les tuples égaux
166 ~~~~~~~~~~~~~~~~~~~~~~~~
167
168 La variable ``d2`` contient 42 tuples mais certains tuples apparaissent plusieurs fois.
169 Pour trouver les tuples égaux à 1e-12 près (précision absolue) invoquer ``findCommonTuples()``.
170 Utiliser ``help(mc.DataArrayDouble.findCommonTuples)`` pour en connaitre l'interface. Stocker le retour de la fonction dans
171 ``c`` et ``cI`` ::
172
173         oldNbOfTuples = d2.getNumberOfTuples()
174         c,cI = d2.findCommonTuples(1e-12)
175
176 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
177 ``i`` dans [0,12). Ainsi les identifiants de tuples du groupe ``i`` commencent à l'index ``cI[i]`` et finissent à l'index ``cI[i+1]``.
178
179 La méthode ``findCommonTuples()`` retourne ainsi 2 paramètres: un tableau contenant la liste des tuples communs 
180 et un tableau d'index qui permet de naviguer dans le premier tableau.    
181 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 
182 maillages non structurés. Cette représentation est rappelée sur l'image ci-dessous, où le premier tableau est en haut, 
183 et le deuxième tableau permettant de la parcourir en bas:
184
185 .. image:: images/IndirectIndex.jpg
186         :scale: 50
187
188
189 .. note:: Le dernier élément de ``cI`` pointe en dehors du tableau ``c``. Ce dernier index est toujours présent 
190    et permet de s'assurer que des traitements tels que les *slices* présentés juste après, sont toujours valables,
191    sans avoir besoin de particulariser le dernier groupe. 
192
193
194 .. _indirect-index-exo:
195
196 Manipuler le format "indirect index"
197 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
198
199 Le nombre de tuples communs à 1e-12 près est donc égal à ``len(cI)-1``, c'est-à-dire 12 dans notre cas.
200 Récupérer la liste des identifiants de tuples du groupe 0 et mettre le résultat dans la variable ``tmp``.
201 Afficher ``tmp``. ::
202
203         tmp = c[cI[0]:cI[0+1]]
204         print tmp
205
206 Vérifier, en l'affichant, que pour tous les identifiants de tuples dans ``tmp``, leurs tuples sont bien égaux dans ``d2``. ::
207
208         print d2[tmp]
209
210 .. 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
211    ces 3 répétitions.
212
213 Maintenant on va déduire des variables ``oldNbOfTuples``, ``c`` et ``cI`` le nombre de tuples effectivement différents dans d2.
214 Pour ce faire, nous allons trouver le nombre de tuples doublons dans ``d2`` et soustraire le résultat de ``oldNbOfTuples``. 
215
216 Pour connaître le nombre de doublons, invoquer ``DataArrayInt.deltaShiftIndex`` qui retourne pour chaque groupe sa taille.
217 Mettre le résultat dans ``a``. ::
218
219         a = cI.deltaShiftIndex()
220
221 Déduire de ``a`` le nombre de tuples doublons dans ``d2`` par groupe et mettre le résultat dans ``b``. ::
222
223         b = a-1
224
225 Enfin on peut trouver le nouveau nombre de tuples grâce à ``b`` et à ``oldNbOfTuples``. Mettre le résultat dans ``myNewNbOfTuples``. ::
226
227         myNewNbOfTuples = oldNbOfTuples - sum(b.getValues())
228
229 Construire un tableau "old-2-new"
230 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231
232 Nous allons maintenant exploiter cette information pour extraire un seul
233 représentant dans chaque groupe de points dupliqués.
234
235 Les deux tableaux ``c`` et ``cI`` définissent une surjection d'un espace de départ à 42 (``oldNbOfTuples``) tuples X 
236 vers un espace à 24 (``myNewNbOfTuples``) tuples Y. 
237
238 .. image:: images/SurjectionDataArray.png
239
240 L'autre manière de définir cette surjection (sans perte d'information) est de la représenter par un tableau "old-2-new".
241 Ce mode de stockage prend la forme d'un DataArrayInt ``o2n`` composé de Card(X) tuples (i.e. 42) à une composante. 
242 Pour chaque tuple (élément) d'index ``i`` de ``o2n``, la case ``o2n[i]`` contient le nouvel identifiant de tuple dans Y.
243 On va donc d'un ancien identifiant (old) vers un nouveau (new).
244
245 Nous allons construire ce tableau pour extraire un sous-ensemble des coordonnées de départ, et ne garder que les 
246 tuples uniques (non doublons) dans l'ensemble de départ.
247
248 .. note:: Pour toutes les opérations de renumérotation en MEDCoupling (bijection), 
249         le format "old-2-new" est systématiquement utilisé.
250
251 La méthode statique ``DataArrayInt.ConvertIndexArrayToO2N()`` (nom un peu barbare, on vous l'accorde) 
252 permet de passer du mode de stockage de cette surjection ``c``, ``cI`` au format ``o2n``.
253 On récupère au passage card(Y) c'est-à-dire le ``newNbOfTuples``. ::
254
255         o2n, newNbOfTuples = mc.DataArrayInt.ConvertIndexArrayToO2N(oldNbOfTuples,c,cI)
256         print "Have I got the right number of tuples?"
257         print "myNewNbOfTuples = %d, newNbOfTuples = %d" % (myNewNbOfTuples, newNbOfTuples)
258         assert(myNewNbOfTuples == newNbOfTuples)
259
260 Nous pouvons maintenant constuire le tableau de points uniques ``d3``. A l'aide de ``o2n`` 
261 et ``newNbOfTuples``, invoquer ``DataArrayDouble.renumberAndReduce()`` sur ``d2``. ::
262
263         d3 = d2.renumberAndReduce(o2n, newNbOfTuples)
264
265 L'inconvénient de cette méthode c'est que finalement on ne connait pas pour chaque groupe de tuple communs dans 
266 d2 quel identifiant a été utilisé.
267 Par exemple pour le groupe 0 on sait que les tuples 0, 8 et 16 (tmp.getValues()) sont tous égaux, et on ne sait 
268 pas si 0, 8 ou 16 a été utilisé pour remplir ``d3``.
269
270 Si l'on souhaite expliciter ce choix, on peut passer en format "new-2-old". Ce mode de stockage prend la forme d'un 
271 ``DataArrayInt`` ``n2o`` composé de Card(Y)
272 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.
273
274 Pour passer d'une description "old-2-new" vers "new-2-old", la méthode est ``DataArrayInt.invertArrayO2N2N2O()``. 
275 Effectuer ce traitement sur la variable ``o2n``. ::
276
277         n2o = o2n.invertArrayO2N2N2O(newNbOfTuples)
278
279 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. ::
280
281         d3_bis = d2[n2o]
282         print "Are d3 and d3_bis equal ? %s" % (str(d3.isEqual(d3_bis, 1e-12)))
283
284 Translater tous les tuples
285 ~~~~~~~~~~~~~~~~~~~~~~~~~~
286
287 Tous les tuples (ou nodes) sont à translater du vecteur [3.3,4.4] afin de recentrer toute la figure en ce point. ::
288
289         d3 += [3.3,4.4]
290
291 Constuire un maillage non structuré
292 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
293
294 On chercher maintenant à créer le maillage final montré dans la figure. Nous avons déjà construit le tableau
295 de coordonnées, il nous reste les cellules à créer.
296
297 Créer un maillage non structuré ``m`` avec les coordonnées ``d3``. Le maillage``m`` a une mesh-dimension 2 ::
298
299         m = mc.MEDCouplingUMesh("My7hexagons",2)
300         m.setCoords(d3)
301         print "Mesh dimension is", m.getMeshDimension()
302         print "Spatial dimension is", m.getCoords().getNumberOfComponents()
303
304 Maintenant, allouer le nombre de cellules avec (un majorant du) nombre attendu de cellules. ::
305
306         m.allocateCells(7)
307
308 Enfin grâce à ``o2n`` on a la *connectivité* (i.e. la liste des points formant un hexagone) 
309 des 7 hexagones utilisant les coordonnées ``d3``. ::
310
311         for i in xrange(7):
312                 cell_connec = o2n[6*i:6*(i+1)]
313                 m.insertNextCell(mc.NORM_POLYGON, cell_connec.getValues())
314                 pass
315
316 Vérifier que ``m`` est correct et ne contient pas d'anomalie. ::
317
318          m.checkConsistencyLight()
319
320 .. note:: Il est toujours une bonne idée d'appeler cette méthode après la construction "from scratch" d'un maillage.
321    Cela assure qu'il n'y a pas de gros "couacs" dans la connectivité, etc ... 
322
323 Pour vérifier *visuellment* que ``m`` est correct, l'écrire dans un fichier "My7hexagons.vtu" et le visualiser dans ParaViS. ::
324
325         m.writeVTK("My7hexagons.vtu")
326
327 .. 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.
328         Bien que l'on écrive au format VTK (\*.vtu), MEDCoupling ne dépend pas de VTK.
329
330 Solution
331 ~~~~~~~~
332
333 :ref:`python_testMEDCouplingdataarray1_solution`