MEDCoupling,  NumPy et SciPy
----------------------------

NumPy est un package additionnel de python qui permet de manipuler des tableaux de manière optimisÊe. 
Il s'agit d'un prĂŠrequis optionnel de MEDCoupling.

NumPy est une passerelle vers le HPC Python (multiprocessing, pyCUDA, SciPy...) et offre de puissantes 
fonctions de calcul vectoriel. C'est pourquoi MEDCoupling offre des liens avec NumPy. 
Un bon point de dĂŠpart pour la dĂŠcouverte de NumPy est le `Tutorial NumPy <http://wiki.scipy.org/Tentative_NumPy_Tutorial>`_

SciPy est aussi un package de python nĂŠcessitant NumPy. Il s'agit ĂŠgalement d'un prĂŠrequis optionnel de MEDCoupling.
SciPy offre des services d'algèbre linÊaire, Fast Fourrier Transform, etc ...

Nous allons ici faire quelques petites manipulations pour voir ce lien entre MEDCoupling et NumPy / SciPy.

DĂŠbut de l'implĂŠmentation
~~~~~~~~~~~~~~~~~~~~~~~~~

Pour commencer l'exercice importer le module Python ``MEDCoupling``: ::

	import MEDCoupling as mc

NumPy est un prĂŠrequis optionnel, vĂŠrifions que nous en bĂŠnĂŠficions bien : ::

	assert(mc.MEDCouplingHasNumPyBindings())

Nous pouvons alors importer NumPy sans problème: ::

	import numpy as np

Convertir un DataArray en tableau NumPy et vice versa
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

CrĂŠer une instance de ``DataArrayDouble`` ayant une composante et 12 tuples.
Et assigner 4. Ă  tous les tuples. ::

	arr = mc.DataArrayDouble(12)
	arr[:] = 4.

CrĂŠons maintenant un tableau NumPy reposant sur les mĂŞmes donnĂŠes que ``arr``. ::

	nparr = arr.toNumPyArray()

Et afficher ``nparr``. ::

	print nparr.__repr__()
	print nparr.tolist()

Mais est ce qu'on ne nous a pas mystifiĂŠ ? ``arr`` et ``nparr`` partagent-ils le mĂŞme bloc mĂŠmoire ?
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. ::

	nparr[::2] = 7.
	print nparr.__repr__()
	print arr.__repr__()

C'est rigolo ! Mais si je dĂŠtruis ``arr`` (le premier Ă  avoir allouĂŠ la mĂŠmoire) est-ce que ``nparr`` est tuĂŠ aussi ? 
Ne risque-t-on pas le SIGSEGV ?
Testons : ::

	del arr
	import gc; gc.collect()     # Make sure the object has been deleted
	print nparr.__repr__()

OK super. Mais inversement puis je faire une instance de ``DataArrayDouble`` avec ``nparr`` ? Oui, en utilisant le constructeur
qui prend un ``nparray`` en entrĂŠe.
Et afficher le contenu.::

	arr2 = mc.DataArrayDouble(nparr)
	print arr2.__repr__()

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.::

	nparr[:] = 5.
	print nparr.__repr__()
	print arr2.__repr__()

Nous en profitons pour montrer un petit service pratique avec NumPy, Ă  savoir, l'ĂŠcriture optimisĂŠe. 
Ecrivons le contenu binaire de ``nparr`` dans un fichier. ::

	f = open("toto.data","w+b")
	a = np.memmap(f, dtype='float64', mode='w+', offset=0, shape=nparr.shape)
	a[:] = nparr[:]
	f.flush()

Relisons "toto.data". ::

	f2 = open("toto.data","r+b")
	b = np.memmap(f2,dtype='float64',mode='r',offset=0,shape=(12,))

Pour rigoler, assignons 3.14 Ă  ``a``, flushons et relisons. ::

	a[:] = 3.14
	f.flush()
	b = np.memmap(f2,dtype='float64',mode='r',offset=0,shape=(12,))
	print b.__repr__()

On voit donc que le passage de MEDCoupling à NumPy se fait directement et de manière optimisÊe. Donc ca peut valoir le coup !
Tout ce qui vient d'ĂŞtre montrĂŠ marche aussi avec des ``DataArrayInt``.
Regardons la suite.

Jouons avec SciPy
~~~~~~~~~~~~~~~~~

Nous allons crĂŠer un maillage non conforme. Le but sera de trouver la peau de ce maillage *sans* les surfaces non conformes.

Nous allons faire cela en jouant avec les matrices creuses de SciPy (*sparse matrix*). Nous interpolons ce maillage non conforme
sur lui mĂŞme, ce qui devrait donner une matrice diagonale si le maillage ĂŠtait conforme.

Avant nous vĂŠrifions que l'on peut jouer avec SciPy ! ::

	assert(mc.MEDCouplingHasSciPyBindings())

Pour le moment crĂŠons un maillage non conforme. Nous collons simplement deux maillages structurĂŠs avec des 
discrĂŠtisations spatiales diffĂŠrentes.::

	c1 = mc.MEDCouplingCMesh()
	arr1 = mc.DataArrayDouble(7) 
	arr1.iota() 
	c1.setCoords(arr1,arr1,arr1)
	c2 = mc.MEDCouplingCMesh()
	arr2 = mc.DataArrayDouble(9)
	arr2.iota() 
	arr2 *= 6./8.
	c2.setCoords(arr2,arr2,arr2)

DĂŠgĂŠnĂŠrons ``c1`` et ``c2`` en non-structurĂŠ, une translation de ``[6.,0.,0.]`` de ``c2``,  et en faisant 
l'agrĂŠgation des deux, c'est dans la poche. ::

	c1 = c1.buildUnstructured()
	c2 = c2.buildUnstructured()
	c2.translate([6.,0.,0.])
	c = mc.MEDCouplingUMesh.MergeUMeshes([c1,c2])

Attention des noeuds sont dupliquĂŠs, il faut invoquer ``mergeNodes()``. ::

	c.mergeNodes(1e-12)

RĂŠcupĂŠrons la peau et les faces non conformes. Ca nous savons faire, car nous avons fait les exercices avant :-) ::

	skinAndNCFaces = c.computeSkin()

Retirons les noeuds non utilisĂŠs. Cette ĂŠtape n'est pas obligatoire. ::

	skinAndNCFaces.zipCoords()

Voici Ă  quoi cela ressemble:

.. image:: images/skinandnccells_numpy.png

OK maintenant on va sÊparer les cellules de bord des cellules non conformes grâce au ``MEDCouplingRemapper``.
Interpolons ``skinAndNCFaces`` sur lui-mĂŞme. On acceptera un ĂŠcart entre face de 1e-12 et un warping max de 0.01. ::

	from MEDCouplingRemapper import MEDCouplingRemapper
	rem = MEDCouplingRemapper()
	rem.setMaxDistance3DSurfIntersect(1e-12)
	rem.setMinDotBtwPlane3DSurfIntersect(0.99)
	rem.prepare(skinAndNCFaces,skinAndNCFaces,"P0P0")

RĂŠcupĂŠrer la matrice creuse au format CSR du remapper. ::

	mat = rem.getCrudeCSRMatrix()
	
.. note:: Le format CSR est un format de stockage efficace des matrices 
	creuses : `Sparse matrix CSR <http://en.wikipedia.org/wiki/Sparse_matrix>`_

Comme nous avons bien suivi les exos sur NumPy, grâce au NumPy array ``mat.indptr`` on peut rÊcupÊrer 
l'ensemble des lignes de la matrice ``mat`` ayant exactement un ĂŠlĂŠment non nul. ::

	indptr = mc.DataArrayInt(mat.indptr)
	indptr2 = indptr.deltaShiftIndex()
	cellIdsOfSkin = indptr2.findIdsEqual(1)

C'est presque fini. CrĂŠer le sous maillage contenant uniquement la peau et l'ĂŠcrire dans 
un fichier VTK ou MED pour le visualiser avec ParaView. ::

	skin = skinAndNCFaces[cellIdsOfSkin]
	skin.writeVTK("skin.vtu")

.. note:: ``skin`` contient des noeuds orphelins, on peut les retirer avec ``skin.zipCoords()``.

Et voilĂ  ce que cela donne :

.. image:: images/skinonly_numpy.png

Script complet
~~~~~~~~~~~~~~

:ref:`python_testMEDCouplingNumPy_solution`