Agitateur - Swirler

import MEDLoader as ml
import numpy as np

# Get available time steps
data = ml.MEDFileData("agitateur.med")
ts = data.getFields()[0].getTimeSteps()
print ts
# Get position of the swirler
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]
# Extract pression field on the swirler
pressMts = data.getFields()["PRESSION_ELEM_DOM"]
press1ts = pressMts[(2,-1)]
pressMc = press1ts.getFieldAtLevel(ml.ON_CELLS,0)
pressOnAgitateurMc = pressMc[ids]
#
pressOnAgitateurMc.getMesh().zipCoords()
# Compute pressure on skin
agitateurMesh3DMc = pressOnAgitateurMc.getMesh()
m3DSurf,desc,descI,revDesc,revDescI = agitateurMesh3DMc.buildDescendingConnectivity()
nbOf3DCellSharing = revDescI.deltaShiftIndex()
ids2 = nbOf3DCellSharing.findIdsEqual(1)
agitateurSkinMc = m3DSurf[ids2]
offsetsOfTupleIdsInField = revDescI[ids2]
tupleIdsInField = revDesc[offsetsOfTupleIdsInField]
pressOnSkinAgitateurMc = pressOnAgitateurMc[tupleIdsInField]
pressOnSkinAgitateurMc.setMesh(agitateurSkinMc)
# Force field computation
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
# Torque computation
singlePolyhedron = agitateurMesh3DMc.buildSpreadZonesWithPoly()
singlePolyhedron.orientCorrectlyPolyhedrons()
centerOfMass = singlePolyhedron.computeCellCenterOfMass()

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

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

zeTorque = torquePerCellOnSkin.accumulate()
print "couple = %r N.m" % zeTorque[2]
# Power computation
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)
# Eigen vector computation
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

def computeAngle(locAgitateur1ts):
        fMc = locAgitateur1ts.getFieldAtLevel(ml.ON_CELLS,0)
        arr = fMc.getArray()
        ids = arr.findIdsInRange(0.,1.)
        f2Mc = fMc[ids]
        m3DSurf,desc,descI,revDesc,revDescI = f2Mc.getMesh().buildDescendingConnectivity()
        nbOf3DCellSharing = revDescI.deltaShiftIndex()
        ids2 = nbOf3DCellSharing.findIdsEqual(1)
        agitateurSkinMc = m3DSurf[ids2]
        #
        singlePolyhedron = agitateurMesh3DMc.buildSpreadZonesWithPoly()
        singlePolyhedron.orientCorrectlyPolyhedrons()
        centerOfMass = singlePolyhedron.computeCellCenterOfMass()
        bary = agitateurSkinMc.computeCellCenterOfMass()
        posSkin = bary-centerOfMass
        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]
        return vect0

vects = len(ts)*[None]
for itts,locAgitateur1ts in zip(ts,data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]):
        angle = computeAngle(locAgitateur1ts)
        vects[itts[0]] = angle
        pass

from math import acos, sqrt
angle2 = len(ts)*[0.]
for pos in xrange(2,len(vects)):
    norm1 = sqrt(vects[pos-1][0]*vects[pos-1][0]+vects[pos-1][1]*vects[pos-1][1])
    norm2 = sqrt(vects[pos][0]*vects[pos][0]+vects[pos][1]*vects[pos][1])
    crs = vects[pos-1][0]*vects[pos][0]+vects[pos-1][1]*vects[pos][1]
    crs /= norm1 ; crs /= norm2 ; crs = min(crs,1.)
    angle2[pos] = acos(crs) #/(ts[pos][2]-ts[pos-1][2])
    pass

omega=sum(angle2)/(ts[-1][2]-ts[0][2])
print sum(angle2)

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)