doc/tutorial/medcouplingloaderex1_fr.rst

   1
   2 MEDCoupling / MEDLoader - Exemple complet 1 - Agitateur
   3 -------------------------------------------------------
   4
   5 Nous partons ici d'un fichier :download:`agitateur.med <data/agitateur.med>` ayant le contenu suivant :
   6
   7 .. image:: images/agitateur.jpg
   8
   9 Il s'agit du résultat d'un petit calcul diphasique : l'agitateur magnétique en vert (repéré seulement par un champ
  10 aux cellules, et n'ayant *pas* de maillage propre) tourne d'un pas de temps à l'autre au
  11 sein d'une phase liquide. Deux gouttes de liquide chutent pendant ce temps vers l'interface air/eau (en gris).
  12
  13 Le but de l'exercice est de calculer le couple appliqué sur cet agitateur, qui est la pièce mécanique entraînant la
  14 partie basse du fluide.
  15
  16 Objectif
  17 ~~~~~~~~
  18
  19 L'objectif est de donner un exemple complet de post-traitement non trivial à partir d'un fichier MED.
  20
  21 Début de l'implémentation
  22 ~~~~~~~~~~~~~~~~~~~~~~~~~
  23
  24 Pour commencer l'exercice importer tout le module python ``MEDLoader`` (qui inclut ``MEDCoupling``).
  25 Importer aussi ``numpy``. ::
  26
  27         import MEDLoader as ml
  28         import numpy as np
  29
  30 Extraction des maillages et champs avec l'API avancée
  31 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  32
  33 Avec l'API avancée lire tout le fichier "agitateur.med" et afficher tous les pas de temps du 1er champ. ::
  34
  35         data = ml.MEDFileData("agitateur.med")
  36         ts = data.getFields()[0].getTimeSteps()
  37         print ts
  38
  39 Récupérer le maillage de l'agitateur (en vert) au pas de temps (2,-1) (cf. ts).
  40 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.
  41 Il faut donc utiliser le champ aux cellules "DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"
  42 et ne sélectionner que la partie du champ ayant une valeur entre dans ``[0.,1.]``. Mettre les identifiants
  43 de cellules correspondant dans ``ids`` : ::
  44
  45         fMts = data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
  46         f1ts = fMts[(2,-1)]
  47         fMc = f1ts.getFieldAtLevel(ml.ON_CELLS,0)
  48         arr = fMc.getArray()
  49         arr.getMinMaxPerComponent()      # just to see the field variation range per component
  50         ids = arr.findIdsInRange(0.,1.)
  51         f2Mc = fMc[ids]
  52
  53 A l'aide du champ "PRESSION_ELEM_DOM" trouver le champ de pression 3D qu'applique l'agitateur. Mettre le résultat dans
  54 ``pressOnAgitateur``. ::
  55
  56         pressMts = data.getFields()["PRESSION_ELEM_DOM"]
  57         press1ts = pressMts[(2,-1)]
  58         pressMc = press1ts.getFieldAtLevel(ml.ON_CELLS,0)
  59         pressOnAgitateurMc = pressMc[ids]
  60
  61 Supprimer les noeuds inutiles de ``pressOnAgitateurMc.getMesh()`` : ::
  62
  63         pressOnAgitateurMc.getMesh().zipCoords()
  64
  65 Passer d'un champ aux cellules 3D à un champ surfacique 3D
  66 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  67
  68 Deduire le champ 3D de pression sur la *peau* de l'agitateur.
  69 Pour ce faire passer par le maillage descendant ``MEDCouplingUMesh.buildDescendingConnectivity()``. ::
  70
  71         agitateurMesh3DMc = pressOnAgitateurMc.getMesh()
  72         m3DSurf,desc,descI,revDesc,revDescI = agitateurMesh3DMc.buildDescendingConnectivity()
  73         nbOf3DCellSharing = revDescI.deltaShiftIndex()
  74         ids2 = nbOf3DCellSharing.findIdsEqual(1)            # Cells with only one neighbor are on the boundary, i.e. on the skin
  75         agitateurSkinMc = m3DSurf[ids2]
  76         offsetsOfTupleIdsInField = revDescI[ids2]
  77         tupleIdsInField = revDesc[offsetsOfTupleIdsInField]
  78         pressOnSkinAgitateurMc = pressOnAgitateurMc[tupleIdsInField]
  79         pressOnSkinAgitateurMc.setMesh(agitateurSkinMc)
  80
  81 Manipuler les champs
  82 ~~~~~~~~~~~~~~~~~~~~
  83
  84 Calculer le champ vectoriel de force sur la peau de l'agitateur en multipliant pour chaque cellule
  85 la pression par la surface et ensuite par le vecteur normal.
  86 La pression est en bar, la convertir au préalable en pascal (Pa). ::
  87
  88         pressSkin = pressOnSkinAgitateurMc.getArray()
  89         pressSkin *= 1e5                   # conversion from bar to Pa
  90         areaSkin = agitateurSkinMc.getMeasureField(True).getArray()
  91         forceSkin = pressSkin*areaSkin
  92         normalSkin = agitateurSkinMc.buildOrthogonalField().getArray()
  93         forceVectSkin = forceSkin*normalSkin
  94
  95 Voici maintenant le premier calcul du moment au centre de masse de l'agitateur :
  96
  97 Pour faire ce 1er calcul de couple exercé sur l'agitateur, calculons la position du centre de masse de l'agitateur.
  98 Calculer le polyèdre représentant l'enveloppe du maillage 3D de l'agitateur ``agitateurMesh3DMc``
  99 (utiliser ``MEDCouplingUMesh.buildSpreadZonesWithPoly()``). ::
 100
 101         singlePolyhedron = agitateurMesh3DMc.buildSpreadZonesWithPoly()
 102         singlePolyhedron.orientCorrectlyPolyhedrons()
 103         centerOfMass = singlePolyhedron.computeCellCenterOfMass()
 104
 105 .. note:: L'appel à ``MEDCouplingUMesh.orientCorrectlyPolyhedrons()`` n'est pas obligatoire mais conseillé car
 106         si par malheur le polyhèdre est mal orienté, son barycentre sera incorrect !
 107
 108 Calculer pour chaque cellule de la peau de l'agitateur le moment par rapport au centre de masse ``centerOfMass``
 109 de l'agitateur.
 110 Pour ce faire calculer ``posSkin`` le ``DataArrayDouble`` donnant pour chaque cellule de la peau de l'agitateur
 111 le vecteur ``centerOfMass`` -> ``G``, avec ``G`` le barycentre de la cellule courante. ::
 112
 113         barySkin=agitateurSkinMc.computeCellCenterOfMass()
 114         posSkin = barySkin-centerOfMass
 115
 116 Appliquer maintenant la formule classique de calcul du moment : calculer le produit
 117 vectoriel par cellule de ``posSkin`` avec ``forceVectSkin`` (méthode ``DataArrayDouble.CrossProduct()``). ::
 118
 119         torquePerCellOnSkin = ml.DataArrayDouble.CrossProduct(posSkin,forceVectSkin)
 120
 121 Sommer ``torqueOnSkin`` en utilisant la méthode ``DataArrayDouble.accumulate()``. ::
 122
 123         zeTorque = torquePerCellOnSkin.accumulate()
 124         print "couple = %r N.m" % zeTorque[2]
 125
 126 Vérifions le couple calculé précédemment en divisant la puissance par la vitesse *angulaire*.
 127 La vitesse *linéaire* est stockée dans le champ "VITESSE_ELEM_DOM".
 128
 129 Calculer la puissance par cellule de la peau de l'agitateur et la sommer. ::
 130
 131         speedMts = data.getFields()["VITESSE_ELEM_DOM"]
 132         speed1ts = speedMts[(2,-1)]
 133         speedMc = speed1ts.getFieldAtLevel(ml.ON_CELLS,0)
 134         speedOnSkin = speedMc.getArray()[tupleIdsInField]
 135         powerSkin = ml.DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
 136         power = powerSkin.accumulate()[0]
 137         print "power = %r W"%(power)
 138
 139 Calculer la vitesse *angulaire*. Pour ce faire, calculer la somme de ``x^2``, ``y^2`` et ``xz`` de ``posSkin`` et
 140 construire (avec NumPy) la matrice 2x2 d'inertie ``inertiaSkin=[[x2,xy], [xy,z2]]``.
 141
 142 Récupérer le vecteur propre associé à la valeur propre maximale
 143 avec ``linalg.eig(inertiaSkin)``. ::
 144
 145         x2 = posSkin[:,0]*posSkin[:,0]
 146         x2 = x2.accumulate()[0]
 147         y2 = posSkin[:,1]*posSkin[:,1]
 148         y2 = y2.accumulate()[0]
 149         xy = posSkin[:,0]*posSkin[:,1]
 150         xy = xy.accumulate()[0]
 151         inertiaSkin = np.matrix([[x2,xy],[xy,y2]])
 152         inertiaSkinValues, inertiaSkinVects = np.linalg.eig(inertiaSkin)
 153         pos = max(enumerate(inertiaSkinValues), key=lambda x: x[1])[0]
 154         vect0 = inertiaSkinVects[pos].tolist()[0]
 155         print vect0
 156
 157 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
 158 détail - on remet les étapes précédentes dans une fonction que l'on applique sur plusieurs pas de temps).
 159
 160 Calculer et comparer le couple sur l'agitateur. ::
 161
 162         omega = 1.1183827931 / (ts[-1][2]-ts[0][2])
 163         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)
 164
 165 Solution
 166 ~~~~~~~~
 167
 168 :ref:`python_testmedcouplingloaderex1_solution`