MEDCoupling / MEDLoader - Exemple complet 1 - Agitateur
-------------------------------------------------------

Nous partons ici d'un fichier :download:`agitateur.med <data/agitateur.med>` ayant le contenu suivant :

.. image:: images/agitateur.jpg

Il s'agit du rĂŠsultat d'un petit calcul diphasique : l'agitateur magnĂŠtique en vert (repĂŠrĂŠ seulement par un champ 
aux cellules, et n'ayant *pas* de maillage propre) tourne d'un pas de temps Ă  l'autre au 
sein d'une phase liquide. Deux gouttes de liquide chutent pendant ce temps vers l'interface air/eau (en gris).  

Le but de l'exercice est de calculer le couple appliquÊ sur cet agitateur, qui est la pièce mÊcanique entraÎnant la
partie basse du fluide.

Objectif
~~~~~~~~

L'objectif est de donner un exemple complet de post-traitement non trivial Ă  partir d'un fichier MED.

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

Pour commencer l'exercice importer tout le module python ``MEDLoader`` (qui inclut ``MEDCoupling``). 
Importer aussi ``numpy``. ::

	import MEDLoader as ml
	import numpy as np

Extraction des maillages et champs avec l'API avancĂŠe
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Avec l'API avancĂŠe lire tout le fichier "agitateur.med" et afficher tous les pas de temps du 1er champ. ::

	data = ml.MEDFileData("agitateur.med")
	ts = data.getFields()[0].getTimeSteps()
	print ts

RĂŠcupĂŠrer le maillage de l'agitateur (en vert) au pas de temps (2,-1) (cf. ts).
La position de l'agitateur est dÊfinie par un champ sur le maillage global du système et n'a pas de maillage propre.
Il faut donc utiliser le champ aux cellules "DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"
et ne sĂŠlectionner que la partie du champ ayant une valeur entre dans ``[0.,1.]``. Mettre les identifiants
de cellules correspondant dans ``ids`` : ::

	fMts = data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
	f1ts = fMts[(2,-1)]
	fMc = f1ts.getFieldAtLevel(ml.ON_CELLS,0)
	arr = fMc.getArray()
	arr.getMinMaxPerComponent()      # just to see the field variation range per component
	ids = arr.findIdsInRange(0.,1.)
	f2Mc = fMc[ids]

A l'aide du champ "PRESSION_ELEM_DOM" trouver le champ de pression 3D qu'applique l'agitateur. Mettre le rĂŠsultat dans
``pressOnAgitateur``. ::

	pressMts = data.getFields()["PRESSION_ELEM_DOM"]
	press1ts = pressMts[(2,-1)]
	pressMc = press1ts.getFieldAtLevel(ml.ON_CELLS,0)
	pressOnAgitateurMc = pressMc[ids]

Supprimer les noeuds inutiles de ``pressOnAgitateurMc.getMesh()`` : ::

	pressOnAgitateurMc.getMesh().zipCoords()

Passer d'un champ aux cellules 3D Ă  un champ surfacique 3D
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Deduire le champ 3D de pression sur la *peau* de l'agitateur.
Pour ce faire passer par le maillage descendant ``MEDCouplingUMesh.buildDescendingConnectivity()``. ::

	agitateurMesh3DMc = pressOnAgitateurMc.getMesh()
	m3DSurf,desc,descI,revDesc,revDescI = agitateurMesh3DMc.buildDescendingConnectivity()
	nbOf3DCellSharing = revDescI.deltaShiftIndex()
	ids2 = nbOf3DCellSharing.findIdsEqual(1)            # Cells with only one neighbor are on the boundary, i.e. on the skin
	agitateurSkinMc = m3DSurf[ids2]
	offsetsOfTupleIdsInField = revDescI[ids2]
	tupleIdsInField = revDesc[offsetsOfTupleIdsInField]
	pressOnSkinAgitateurMc = pressOnAgitateurMc[tupleIdsInField]
	pressOnSkinAgitateurMc.setMesh(agitateurSkinMc)

Manipuler les champs
~~~~~~~~~~~~~~~~~~~~

Calculer le champ vectoriel de force sur la peau de l'agitateur en multipliant pour chaque cellule
la pression par la surface et ensuite par le vecteur normal.
La pression est en bar, la convertir au prĂŠalable en pascal (Pa). ::

	pressSkin = pressOnSkinAgitateurMc.getArray()
	pressSkin *= 1e5                   # conversion from bar to Pa
	areaSkin = agitateurSkinMc.getMeasureField(True).getArray()  
	forceSkin = pressSkin*areaSkin
	normalSkin = agitateurSkinMc.buildOrthogonalField().getArray()
	forceVectSkin = forceSkin*normalSkin

Voici maintenant le premier calcul du moment au centre de masse de l'agitateur :

Pour faire ce 1er calcul de couple exercĂŠ sur l'agitateur, calculons la position du centre de masse de l'agitateur.
Calculer le polyèdre reprÊsentant l'enveloppe du maillage 3D de l'agitateur ``agitateurMesh3DMc``
(utiliser ``MEDCouplingUMesh.buildSpreadZonesWithPoly()``). ::

	singlePolyhedron = agitateurMesh3DMc.buildSpreadZonesWithPoly()
	singlePolyhedron.orientCorrectlyPolyhedrons()
	centerOfMass = singlePolyhedron.computeCellCenterOfMass()

.. note:: L'appel Ă  ``MEDCouplingUMesh.orientCorrectlyPolyhedrons()`` n'est pas obligatoire mais conseillĂŠ car 
	si par malheur le polyhèdre est mal orientÊ, son barycentre sera incorrect !

Calculer pour chaque cellule de la peau de l'agitateur le moment par rapport au centre de masse ``centerOfMass``
de l'agitateur.
Pour ce faire calculer ``posSkin`` le ``DataArrayDouble`` donnant pour chaque cellule de la peau de l'agitateur
le vecteur ``centerOfMass`` -> ``G``, avec ``G`` le barycentre de la cellule courante. ::

	barySkin=agitateurSkinMc.computeCellCenterOfMass()
	posSkin = barySkin-centerOfMass

Appliquer maintenant la formule classique de calcul du moment : calculer le produit 
vectoriel par cellule de ``posSkin`` avec ``forceVectSkin`` (mĂŠthode ``DataArrayDouble.CrossProduct()``). ::

	torquePerCellOnSkin = ml.DataArrayDouble.CrossProduct(posSkin,forceVectSkin)

Sommer ``torqueOnSkin`` en utilisant la mĂŠthode ``DataArrayDouble.accumulate()``. ::

	zeTorque = torquePerCellOnSkin.accumulate()
	print "couple = %r N.m" % zeTorque[2]

VĂŠrifions le couple calculĂŠ prĂŠcĂŠdemment en divisant la puissance par la vitesse *angulaire*.
La vitesse *linĂŠaire* est stockĂŠe dans le champ "VITESSE_ELEM_DOM".

Calculer la puissance par cellule de la peau de l'agitateur et la sommer. ::

	speedMts = data.getFields()["VITESSE_ELEM_DOM"]
	speed1ts = speedMts[(2,-1)]
	speedMc = speed1ts.getFieldAtLevel(ml.ON_CELLS,0)
	speedOnSkin = speedMc.getArray()[tupleIdsInField]
	powerSkin = ml.DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
	power = powerSkin.accumulate()[0]
	print "power = %r W"%(power)

Calculer la vitesse *angulaire*. Pour ce faire, calculer la somme de ``x^2``, ``y^2`` et ``xz`` de ``posSkin`` et 
construire (avec NumPy) la matrice 2x2 d'inertie ``inertiaSkin=[[x2,xy], [xy,z2]]``.

RĂŠcupĂŠrer le vecteur propre associĂŠ Ă  la valeur propre maximale
avec ``linalg.eig(inertiaSkin)``. ::

	x2 = posSkin[:,0]*posSkin[:,0]
	x2 = x2.accumulate()[0]
	y2 = posSkin[:,1]*posSkin[:,1]
	y2 = y2.accumulate()[0]
	xy = posSkin[:,0]*posSkin[:,1]
	xy = xy.accumulate()[0]
	inertiaSkin = np.matrix([[x2,xy],[xy,y2]])
	inertiaSkinValues, inertiaSkinVects = np.linalg.eig(inertiaSkin)
	pos = max(enumerate(inertiaSkinValues), key=lambda x: x[1])[0]
	vect0 = inertiaSkinVects[pos].tolist()[0]
	print vect0

Grâce au calcul prÊcÊdent on peut dÊduire que l'agitateur a tournÊ de 1.1183827931 radian (cf. solution complète pour le
dĂŠtail - on remet les ĂŠtapes prĂŠcĂŠdentes dans une fonction que l'on applique sur plusieurs pas de temps).

Calculer et comparer le couple sur l'agitateur. ::

	omega = 1.1183827931 / (ts[-1][2]-ts[0][2])
	print "At timestep (%d,%d) (physical time=%r s) the torque is: %r N.m, power/omega=%r N.m " % (ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega)

Solution
~~~~~~~~

:ref:`python_testmedcouplingloaderex1_solution`