MEDCouplingRemapper : interpolation de champs
---------------------------------------------

Nous allons ici effectuer une interpolation entre deux maillages ``srcMesh`` et ``trgMesh``.
Pour mettre l'accent sur certaines petites subtilitĂŠs de l'interpolation, nous prenons un cas particulier oĂš ``srcMesh`` est 
un maillage raffinĂŠ de ``trgMesh`` (avec certaines cellules dĂŠcoupĂŠes plus finement).

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

Pour commencer l'exercice importer le module Python ``MEDCoupling`` et la 
classe ``MEDCouplingRemapper`` du module ``MEDCouplingRemapper``. ::

	import MEDCoupling as mc
	from MEDCouplingRemapper import MEDCouplingRemapper 


CrĂŠer le maillage target
~~~~~~~~~~~~~~~~~~~~~~~~

Construire le maillage non structurÊ ``trgMesh`` issu d'un maillage cartÊsien 2D 10x10 commençant au 
point ``[0.,0.]`` et ayant un pas de 1. selon X et selon Y : ::

	arr = mc.DataArrayDouble(11)
	arr.iota(0)
	trgMesh = mc.MEDCouplingCMesh()
	trgMesh.setCoords(arr,arr)
	trgMesh = trgMesh.buildUnstructured()	

CrĂŠer le maillage source
~~~~~~~~~~~~~~~~~~~~~~~~

CrÊer un maillage ``srcMesh`` issu d'un maillage cartÊsien 2D 20x20 cellules commençant 
aussi au point ``[0.,0.]`` et ayant un pas de 0.5 selon X et selon Y : ::

	arr = mc.DataArrayDouble(21)
	arr.iota(0)
	arr *= 0.5
	srcMesh = mc.MEDCouplingCMesh()
	srcMesh.setCoords(arr,arr)
	srcMesh = srcMesh.buildUnstructured()	

Afin de rendre l'exercise plus intĂŠressant, triangulariser Ă  l'aide de ``MEDCouplingUMesh.simplexize()`` 
les 20 premières cellules de ``srcMesh``
(les simplexes 2D sont les triangles). Mettre le rĂŠsultat dans ``srcMesh``. ::

	tmp = srcMesh[:20]    # Extract a sub-part of srcMesh
	tmp.simplexize(0)
	srcMesh = mc.MEDCouplingUMesh.MergeUMeshes([tmp,srcMesh[20:]])

Interpoler avec MEDCouplingRemapper
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Nous rappelons que pour projeter un champ d'un maillage vers un autre, il faut d'abord prĂŠparer la matrice d'interpolation
qui contient les ratios de projection.

Calculer la première partie de la matrice d'interpolation de ``srcMesh`` (discrÊtisÊe aux cellules - *P0*) 
vers ``trgMesh`` (discrĂŠtisĂŠe aux cellules ĂŠgalement).
Pour ce faire, invoquer ``MEDCouplingRemapper.prepare()`` sur une instance (``remap``) de ``MEDCouplingRemapper``. ::

	remap = MEDCouplingRemapper()
	remap.prepare(srcMesh,trgMesh,"P0P0")

VĂŠrifier que la matrice calculĂŠe par la mĂŠthode est correcte dans notre cas trivial.
Pour ce faire, rĂŠcupĂŠrer dans ``myMatrix`` la matrice interne retournĂŠe par ``MEDCouplingRemapper.getCrudeMatrix()``.
Celle-ci donne pour chaque cellule de ``trgMesh`` les identifiants de cellules de ``srcMesh`` avec
lesquelles elle s'intersecte, et l'aire d'intersection correspondante. 

VĂŠrifier notamment que pour chaque cellule de ``trgMesh`` la somme des aires fait toujours 1. ::

	myMatrix = remap.getCrudeMatrix()
	print myMatrix
	sumByRows = mc.DataArrayDouble(len(myMatrix))
	for i,wIt in enumerate(sumByRows):
	  su = 0.
	  for it in myMatrix[i]:
	    su += myMatrix[i][it]
	  wIt[0] = su
	print "Is interpolation well prepared?", sumByRows.isUniform(1.,1e-12)

.. note:: Les triangles dans ``srcMesh`` ont ĂŠtĂŠ rajoutĂŠs pour casser la monotonie de la matrice ``myMatrix``.

.. note:: Comme on le voit, la prĂŠparation ne nĂŠcessite que les deux maillages, et rien d'autre.

Construire un champ aux cellules "srcField" construit Ă  partir de la formule analytique 
suivante ``7-sqrt((x-5.)*(x-5.)+(y-5.)*(y-5.))`` : ::

	srcField = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
	srcField.setMesh(srcMesh)
	srcField.fillFromAnalytic(1,"7-sqrt((x-5.)*(x-5.)+(y-5.)*(y-5.))")
	srcField.getArray().setInfoOnComponent(0, "powercell [W]")

Voici Ă  quoi ressemble ce champ :

.. image:: images/Remapper1.png

Appliquer l'interpolation avec ``MEDCouplingRemapper.transferField()`` : ::

	remap.transferField(srcField, 1e300)

.. note:: 1e300 est une valeur par dĂŠfaut. Cette valeur sera systĂŠmatiquement assignĂŠe Ă  toute cellule
	de ``trgField`` n'interceptant aucune cellule de ``srcMesh``. En gĂŠnĂŠral, les utilisateurs mettent une 
	valeur ĂŠnorme pour repĂŠrer ce qui est souvent un bug. Mais d'autres utilisateurs, dans la perspective 
	d'une interpolation parallèle par exemple, mettent 0.

.. note:: Une exception est envoyĂŠe car ``srcField`` n'a pas de *nature* dĂŠfinie. 
	Nous allons voir dans la suite l'impact de cet attribut sur le rĂŠsultat final.

Mettre la nature de ``srcField`` Ă  ``IntensiveMaximum``. Cela signifie que le champ doit ĂŞtre interprĂŠtĂŠ commĂŠ ĂŠtant
intensif (une tempĂŠrature par exemple). ::

	srcField.setNature(mc.IntensiveMaximum)
	trgFieldCV = remap.transferField(srcField,1e300)

VĂŠrifier qu'avec la nature ``IntensiveMaximum``, l'intĂŠgrale du champ est conservĂŠe. Par contre, 
la somme sur les cellules (accumulation) n'est **pas** conservĂŠe ! ::

	integSource = srcField.integral(True)[0]
	integTarget =  trgFieldCV.integral(True)[0]
	print "IntensiveMaximum -- integrals: %lf == %lf" % (integSource, integTarget)
	
	accSource = srcField.getArray().accumulate()[0]
	accTarget = trgFieldCV.getArray().accumulate()[0]
	print "IntensiveMaximum -- sums: %lf != %lf" % (accSource, accTarget)


Maintenant mettre la nature de ``srcField`` Ă  ``ExtensiveConservation``. Le champ doit ĂŞtre interprĂŠtĂŠ commĂŠ ĂŠtant
extensif (par exemple une puissance ou un volume). ::

	srcField.setNature(mc.ExtensiveConservation)
	trgFieldI = remap.transferField(srcField,1e300)

VĂŠrifier qu'avec la nature ``ExtensiveConservation``, l'intĂŠgrale du champ n'est **pas** conservĂŠe. 
Par contre, la somme sur les cellules est conservĂŠe. ::

	integSource = srcField.integral(True)[0]
	integTarget =  trgFieldI.integral(True)[0]
	print "ExtensiveConservation -- integrals: %lf != %lf" % (integSource, integTarget)
	
	accSource = srcField.getArray().accumulate()[0]
	accTarget = trgFieldI.getArray().accumulate()[0]
	print "ExtensiveConservation -- sums: %lf == %lf" % (accSource, accTarget)

Visualiser les champs avec ParaViS, ou en les ĂŠcrivant dans un fichier.

Solution
~~~~~~~~

:ref:`python_testMEDCouplingremapper1_solution`