]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/atestMEDCouplingLoaderEx1.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / 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)