]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/medcouplingnumpy_en.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / doc / tutorial / medcouplingnumpy_en.rst
1
2 MEDCoupling,  NumPy et SciPy
3 ----------------------------
4
5 NumPy est un package additionnel de python qui permet de manipuler des tableaux de manière optimisée. 
6 Il s'agit d'un prérequis optionnel de MEDCoupling.
7
8 NumPy est une passerelle vers le HPC Python (multiprocessing, pyCUDA, SciPy...) et offre de puissantes 
9 fonctions de calcul vectoriel. C'est pourquoi MEDCoupling offre des liens avec NumPy. 
10 Un bon point de départ pour la découverte de NumPy est le `Tutorial NumPy <http://wiki.scipy.org/Tentative_NumPy_Tutorial>`_
11
12 SciPy est aussi un package de python nécessitant NumPy. Il s'agit également d'un prérequis optionnel de MEDCoupling.
13 SciPy offre des services d'algèbre linéaire, Fast Fourrier Transform, etc ...
14
15 Nous allons ici faire quelques petites manipulations pour voir ce lien entre MEDCoupling et NumPy / SciPy.
16
17 Début de l'implémentation
18 ~~~~~~~~~~~~~~~~~~~~~~~~~
19
20 Pour commencer l'exercice importer le module Python ``MEDCoupling``: ::
21
22         import MEDCoupling as mc
23
24 NumPy est un prérequis optionnel, vérifions que nous en bénéficions bien : ::
25
26         assert(mc.MEDCouplingHasNumPyBindings())
27
28 Nous pouvons alors importer NumPy sans problème: ::
29
30         import numpy as np
31
32 Convertir un DataArray en tableau NumPy et vice versa
33 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
34
35 Créer une instance de ``DataArrayDouble`` ayant une composante et 12 tuples.
36 Et assigner 4. à tous les tuples. ::
37
38         arr = mc.DataArrayDouble(12)
39         arr[:] = 4.
40
41 Créons maintenant un tableau NumPy reposant sur les mêmes données que ``arr``. ::
42
43         nparr = arr.toNumPyArray()
44
45 Et afficher ``nparr``. ::
46
47         print nparr.__repr__()
48         print nparr.tolist()
49
50 Mais est ce qu'on ne nous a pas mystifié ? ``arr`` et ``nparr`` partagent-ils le même bloc mémoire ?
51 Pour le vérifier assignons 7.0 un tuple sur 2 avec ``nparr`` et vérifions que ``arr`` et ``nparr`` sont simultanément modifiés. ::
52
53         nparr[::2] = 7.
54         print nparr.__repr__()
55         print arr.__repr__()
56
57 C'est rigolo ! Mais si je détruis ``arr`` (le premier à avoir alloué la mémoire) est-ce que ``nparr`` est tué aussi ? 
58 Ne risque-t-on pas le SIGSEGV ?
59 Testons : ::
60
61         del arr
62         import gc; gc.collect()     # Make sure the object has been deleted
63         print nparr.__repr__()
64
65 OK super. Mais inversement puis je faire une instance de ``DataArrayDouble`` avec ``nparr`` ? Oui, en utilisant le constructeur
66 qui prend un ``nparray`` en entrée.
67 Et afficher le contenu.::
68
69         arr2 = mc.DataArrayDouble(nparr)
70         print arr2.__repr__()
71
72 Modifions ``nparr`` en assignant 5.0 pour tous les tuples et vérifier que les 2 représentations ont bien été modifiées simultanément.::
73
74         nparr[:] = 5.
75         print nparr.__repr__()
76         print arr2.__repr__()
77
78 Nous en profitons pour montrer un petit service pratique avec NumPy, à savoir, l'écriture optimisée. 
79 Ecrivons le contenu binaire de ``nparr`` dans un fichier. ::
80
81         f = open("toto.data","w+b")
82         a = np.memmap(f, dtype='float64', mode='w+', offset=0, shape=nparr.shape)
83         a[:] = nparr[:]
84         f.flush()
85
86 Relisons "toto.data". ::
87
88         f2 = open("toto.data","r+b")
89         b = np.memmap(f2,dtype='float64',mode='r',offset=0,shape=(12,))
90
91 Pour rigoler, assignons 3.14 à ``a``, flushons et relisons. ::
92
93         a[:] = 3.14
94         f.flush()
95         b = np.memmap(f2,dtype='float64',mode='r',offset=0,shape=(12,))
96         print b.__repr__()
97
98 On voit donc que le passage de MEDCoupling à NumPy se fait directement et de manière optimisée. Donc ca peut valoir le coup !
99 Tout ce qui vient d'être montré marche aussi avec des ``DataArrayInt``.
100 Regardons la suite.
101
102 Jouons avec SciPy
103 ~~~~~~~~~~~~~~~~~
104
105 Nous allons créer un maillage non conforme. Le but sera de trouver la peau de ce maillage *sans* les surfaces non conformes.
106
107 Nous allons faire cela en jouant avec les matrices creuses de SciPy (*sparse matrix*). Nous interpolons ce maillage non conforme
108 sur lui même, ce qui devrait donner une matrice diagonale si le maillage était conforme.
109
110 Avant nous vérifions que l'on peut jouer avec SciPy ! ::
111
112         assert(mc.MEDCouplingHasSciPyBindings())
113
114 Pour le moment créons un maillage non conforme. Nous collons simplement deux maillages structurés avec des 
115 discrétisations spatiales différentes.::
116
117         c1 = mc.MEDCouplingCMesh()
118         arr1 = mc.DataArrayDouble(7) 
119         arr1.iota() 
120         c1.setCoords(arr1,arr1,arr1)
121         c2 = mc.MEDCouplingCMesh()
122         arr2 = mc.DataArrayDouble(9)
123         arr2.iota() 
124         arr2 *= 6./8.
125         c2.setCoords(arr2,arr2,arr2)
126
127 Dégénérons ``c1`` et ``c2`` en non-structuré, une translation de ``[6.,0.,0.]`` de ``c2``,  et en faisant 
128 l'agrégation des deux, c'est dans la poche. ::
129
130         c1 = c1.buildUnstructured()
131         c2 = c2.buildUnstructured()
132         c2.translate([6.,0.,0.])
133         c = mc.MEDCouplingUMesh.MergeUMeshes([c1,c2])
134
135 Attention des noeuds sont dupliqués, il faut invoquer ``mergeNodes()``. ::
136
137         c.mergeNodes(1e-12)
138
139 Récupérons la peau et les faces non conformes. Ca nous savons faire, car nous avons fait les exercices avant :-) ::
140
141         skinAndNCFaces = c.computeSkin()
142
143 Retirons les noeuds non utilisés. Cette étape n'est pas obligatoire. ::
144
145         skinAndNCFaces.zipCoords()
146
147 Voici à quoi cela ressemble:
148
149 .. image:: images/skinandnccells_numpy.png
150
151 OK maintenant on va séparer les cellules de bord des cellules non conformes grâce au ``MEDCouplingRemapper``.
152 Interpolons ``skinAndNCFaces`` sur lui-même. On acceptera un écart entre face de 1e-12 et un warping max de 0.01. ::
153
154         from MEDCouplingRemapper import MEDCouplingRemapper
155         rem = MEDCouplingRemapper()
156         rem.setMaxDistance3DSurfIntersect(1e-12)
157         rem.setMinDotBtwPlane3DSurfIntersect(0.99)
158         rem.prepare(skinAndNCFaces,skinAndNCFaces,"P0P0")
159
160 Récupérer la matrice creuse au format CSR du remapper. ::
161
162         mat = rem.getCrudeCSRMatrix()
163         
164 .. note:: Le format CSR est un format de stockage efficace des matrices 
165         creuses : `Sparse matrix CSR <http://en.wikipedia.org/wiki/Sparse_matrix>`_
166
167 Comme nous avons bien suivi les exos sur NumPy, grâce au NumPy array ``mat.indptr`` on peut récupérer 
168 l'ensemble des lignes de la matrice ``mat`` ayant exactement un élément non nul. ::
169
170         indptr = mc.DataArrayInt(mat.indptr)
171         indptr2 = indptr.deltaShiftIndex()
172         cellIdsOfSkin = indptr2.findIdsEqual(1)
173
174 C'est presque fini. Créer le sous maillage contenant uniquement la peau et l'écrire dans 
175 un fichier VTK ou MED pour le visualiser avec ParaView. ::
176
177         skin = skinAndNCFaces[cellIdsOfSkin]
178         skin.writeVTK("skin.vtu")
179
180 .. note:: ``skin`` contient des noeuds orphelins, on peut les retirer avec ``skin.zipCoords()``.
181
182 Et voilà ce que cela donne :
183
184 .. image:: images/skinonly_numpy.png
185
186 Script complet
187 ~~~~~~~~~~~~~~
188
189 :ref:`python_testMEDCouplingNumPy_solution`
190