doc/tutorial/atestMEDCouplingLoaderEx1.rst

   1
   2 .. _python_testmedcouplingloaderex1_solution:
   3
   4 Agitateur - Swirler
   5 ~~~~~~~~~~~~~~~~~~~
   6
   7 ::
   8
   9         import MEDLoader as ml
  10         import numpy as np
  11
  12         # Get available time steps
  13         data = ml.MEDFileData("agitateur.med")
  14         ts = data.getFields()[0].getTimeSteps()
  15         print ts
  16         # Get position of the swirler
  17         fMts = data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
  18         f1ts = fMts[(2,-1)]
  19         fMc = f1ts.getFieldAtLevel(ml.ON_CELLS,0)
  20         arr = fMc.getArray()
  21         arr.getMinMaxPerComponent()      # just to see the field variation range per component
  22         ids = arr.findIdsInRange(0.,1.)
  23         f2Mc = fMc[ids]
  24         # Extract pression field on the swirler
  25         pressMts = data.getFields()["PRESSION_ELEM_DOM"]
  26         press1ts = pressMts[(2,-1)]
  27         pressMc = press1ts.getFieldAtLevel(ml.ON_CELLS,0)
  28         pressOnAgitateurMc = pressMc[ids]
  29         #
  30         pressOnAgitateurMc.getMesh().zipCoords()
  31         # Compute pressure on skin
  32         agitateurMesh3DMc = pressOnAgitateurMc.getMesh()
  33         m3DSurf,desc,descI,revDesc,revDescI = agitateurMesh3DMc.buildDescendingConnectivity()
  34         nbOf3DCellSharing = revDescI.deltaShiftIndex()
  35         ids2 = nbOf3DCellSharing.findIdsEqual(1)
  36         agitateurSkinMc = m3DSurf[ids2]
  37         offsetsOfTupleIdsInField = revDescI[ids2]
  38         tupleIdsInField = revDesc[offsetsOfTupleIdsInField]
  39         pressOnSkinAgitateurMc = pressOnAgitateurMc[tupleIdsInField]
  40         pressOnSkinAgitateurMc.setMesh(agitateurSkinMc)
  41         # Force field computation
  42         pressSkin = pressOnSkinAgitateurMc.getArray()
  43         pressSkin *= 1e5                   # conversion from bar to Pa
  44         areaSkin = agitateurSkinMc.getMeasureField(True).getArray()
  45         forceSkin = pressSkin*areaSkin
  46         normalSkin = agitateurSkinMc.buildOrthogonalField().getArray()
  47         forceVectSkin = forceSkin*normalSkin
  48         # Torque computation
  49         singlePolyhedron = agitateurMesh3DMc.buildSpreadZonesWithPoly()
  50         singlePolyhedron.orientCorrectlyPolyhedrons()
  51         centerOfMass = singlePolyhedron.computeCellCenterOfMass()
  52
  53         barySkin=agitateurSkinMc.computeCellCenterOfMass()
  54         posSkin = barySkin-centerOfMass
  55
  56         torquePerCellOnSkin = ml.DataArrayDouble.CrossProduct(posSkin,forceVectSkin)
  57
  58         zeTorque = torquePerCellOnSkin.accumulate()
  59         print "couple = %r N.m" % zeTorque[2]
  60         # Power computation
  61         speedMts = data.getFields()["VITESSE_ELEM_DOM"]
  62         speed1ts = speedMts[(2,-1)]
  63         speedMc = speed1ts.getFieldAtLevel(ml.ON_CELLS,0)
  64         speedOnSkin = speedMc.getArray()[tupleIdsInField]
  65         powerSkin = ml.DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
  66         power = powerSkin.accumulate()[0]
  67         print "power = %r W"%(power)
  68         # Eigen vector computation
  69         x2 = posSkin[:,0]*posSkin[:,0]
  70         x2 = x2.accumulate()[0]
  71         y2 = posSkin[:,1]*posSkin[:,1]
  72         y2 = y2.accumulate()[0]
  73         xy = posSkin[:,0]*posSkin[:,1]
  74         xy = xy.accumulate()[0]
  75         inertiaSkin = np.matrix([[x2,xy],[xy,y2]])
  76         inertiaSkinValues, inertiaSkinVects = np.linalg.eig(inertiaSkin)
  77         pos = max(enumerate(inertiaSkinValues), key=lambda x: x[1])[0]
  78         vect0 = inertiaSkinVects[pos].tolist()[0]
  79         print vect0
  80
  81         def computeAngle(locAgitateur1ts):
  82                 fMc = locAgitateur1ts.getFieldAtLevel(ml.ON_CELLS,0)
  83                 arr = fMc.getArray()
  84                 ids = arr.findIdsInRange(0.,1.)
  85                 f2Mc = fMc[ids]
  86                 m3DSurf,desc,descI,revDesc,revDescI = f2Mc.getMesh().buildDescendingConnectivity()
  87                 nbOf3DCellSharing = revDescI.deltaShiftIndex()
  88                 ids2 = nbOf3DCellSharing.findIdsEqual(1)
  89                 agitateurSkinMc = m3DSurf[ids2]
  90                 #
  91                 singlePolyhedron = agitateurMesh3DMc.buildSpreadZonesWithPoly()
  92                 singlePolyhedron.orientCorrectlyPolyhedrons()
  93                 centerOfMass = singlePolyhedron.computeCellCenterOfMass()
  94                 bary = agitateurSkinMc.computeCellCenterOfMass()
  95                 posSkin = bary-centerOfMass
  96                 x2=posSkin[:,0]*posSkin[:,0] ; x2=x2.accumulate()[0]
  97                 y2=posSkin[:,1]*posSkin[:,1] ; y2=y2.accumulate()[0]
  98                 xy=posSkin[:,0]*posSkin[:,1] ; xy=xy.accumulate()[0]
  99                 inertiaSkin = np.matrix([[x2,xy],[xy,y2]])
 100                 inertiaSkinValues,inertiaSkinVects = np.linalg.eig(inertiaSkin)
 101                 pos = max(enumerate(inertiaSkinValues), key=lambda x: x[1])[0]
 102                 vect0 = inertiaSkinVects[pos].tolist()[0]
 103                 return vect0
 104
 105         vects = len(ts)*[None]
 106         for itts,locAgitateur1ts in zip(ts,data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]):
 107                 angle = computeAngle(locAgitateur1ts)
 108                 vects[itts[0]] = angle
 109                 pass
 110
 111         from math import acos, sqrt
 112         angle2 = len(ts)*[0.]
 113         for pos in xrange(2,len(vects)):
 114             norm1 = sqrt(vects[pos-1][0]*vects[pos-1][0]+vects[pos-1][1]*vects[pos-1][1])
 115             norm2 = sqrt(vects[pos][0]*vects[pos][0]+vects[pos][1]*vects[pos][1])
 116             crs = vects[pos-1][0]*vects[pos][0]+vects[pos-1][1]*vects[pos][1]
 117             crs /= norm1 ; crs /= norm2 ; crs = min(crs,1.)
 118             angle2[pos] = acos(crs) #/(ts[pos][2]-ts[pos-1][2])
 119             pass
 120
 121         omega=sum(angle2)/(ts[-1][2]-ts[0][2])
 122         print sum(angle2)
 123
 124         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)