medtool/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.getIdsEqual(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