doc/tutorial/medcouplingnumpy_fr.rst

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