3 Full example 1 - Agitator
4 -------------------------
6 The MED file :download:`agitateur.med <data/agitateur.med>` is used and has the following content:ant :
8 .. image:: images/agitateur.jpg
10 This is the result of a simple 2-phase computation.
11 The agitator in green (represented by a cell field) turn on itself from one time-step
13 The purpose of the exercise is to compute the torque applied on this piece.agitateur.
18 The aim of this exercise is to give a full example of non-trivial post-treatment
24 Import the whole Python module medcoupling.
27 import medcoupling as mc
30 Mesh and field extraction using advanced API
31 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
33 Using the advanced API read the whole file "agitateur.med" and display all time-steps of
36 data=mc.MEDFileData("agitateur.med")
37 ts=data.getFields()[0].getTimeSteps()
40 Get the agitator's mesh (in green) at the time-step (2,-1) (see ts).
41 To this end use the cell field "DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM" and select
42 only the field part having a value within [0.0, 1.0] (variable "ids"). ::
44 fMts=data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
46 fMc=f1ts.getFieldAtLevel(mc.ON_CELLS,0)
48 arr.getMinMaxPerComponent() # just to see the variation range of the field per component
49 ids=arr.findIdsInRange(0.,1.)
52 Using the field "PRESSION_ELEM_DOM" find the 3D pression field applied on the agitator.
53 Store the result in pressOnAgitateur. ::
55 pressMts=data.getFields()["PRESSION_ELEM_DOM"]
56 press1ts=pressMts[(2,-1)]
57 pressMc=press1ts.getFieldAtLevel(mc.ON_CELLS,0)
58 pressOnAgitateurMc=pressMc[ids]
60 Delete unused nodes in pressOnAgitateurMc.getMesh(). ::
62 pressOnAgitateurMc.getMesh().zipCoords()
64 Create a 3D surface field from the 3D cell field
65 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
67 Deduce the 3D field on the skin of the agitator.
68 To achieve this use the constituting mesh MEDCouplingUMesh.buildDescendingConnectivity().
69 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
71 agitateurMesh3DMc=pressOnAgitateurMc.getMesh()
72 m3DSurf,desc,descI,revDesc,revDescI=agitateurMesh3DMc.buildDescendingConnectivity()
73 nbOf3DCellSharing=revDescI.deltaShiftIndex()
74 ids2=nbOf3DCellSharing.findIdsEqual(1)
75 agitateurSkinMc=m3DSurf[ids2]
76 OffsetsOfTupleIdsInField=revDescI[ids2]
77 tupleIdsInField=revDesc[OffsetsOfTupleIdsInField]
78 pressOnSkinAgitateurMc=pressOnAgitateurMc[tupleIdsInField]
79 pressOnSkinAgitateurMc.setMesh(agitateurSkinMc)
84 Compute the force vector field on the agitator's skin by multiplying at each cell
85 the pressure by the surface, and then the normal vector.
86 Pression is expressed in bar, convert it first to Pa. ::
88 pressSkin=pressOnSkinAgitateurMc.getArray()
90 areaSkin=agitateurSkinMc.getMeasureField(True).getArray()
91 forceSkin=pressSkin*areaSkin
92 normalSkin=agitateurSkinMc.buildOrthogonalField().getArray()
93 forceVectSkin=forceSkin*normalSkin
95 First computation of the torque at the center of mass of the agitator:
97 Let's compute first the position of the center of mass.
98 Compute the polyhedron representing the 3D mesh hull of the agitator "agitateurMesh3DMc"
99 (use MEDCouplingUMesh.buildSpreadZonesWithPoly()). ::
101 singlePolyhedron=agitateurMesh3DMc.buildSpreadZonesWithPoly()
102 singlePolyhedron.orientCorrectlyPolyhedrons()
103 centerOfMass=singlePolyhedron.computeCellCenterOfMass()
105 .. note:: The call to MEDCouplingUMesh.orientCorrectlyPolyhedrons() is not mandatory
106 but is recommended: if the polyhedron happens to be mis-oriented, its center of mass will
109 Compute for each skin cell the torque with respect to the center of mass "centerOfMass".
110 To this end compute "posSkin", a DataArrayDouble giving for each skin cell the vector
111 centerOfMass -> G, where G represents the center of mass of the current cell. ::
113 barySkin=agitateurSkinMc.computeCellCenterOfMass()
114 posSkin=barySkin-centerOfMass
116 Compute the cross product for each cell of "posSkin" using "forceVectSkin"
117 (method DataArrayDouble.CrossProduct()). ::
119 torquePerCellOnSkin=mc.DataArrayDouble.CrossProduct(posSkin,forceVectSkin)
121 Sum "torqueOnSkin" using DataArrayDouble.accumulate(). ::
123 zeTorque=torquePerCellOnSkin.accumulate()
124 print("couple = %r N.m"%(zeTorque[2]))
126 Check the previously computed torque by dividing the power by the angular speed.
127 Compute the power per skin cell and sum it. ::
129 speedMts=data.getFields()["VITESSE_ELEM_DOM"]
130 speed1ts=speedMts[(2,-1)]
131 speedMc=speed1ts.getFieldAtLevel(mc.ON_CELLS,0)
132 speedOnSkin=speedMc.getArray()[tupleIdsInField]
133 powerSkin=mc.DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
134 power=powerSkin.accumulate()[0]
135 print("power = %r W"%(power))
137 Compute the angular speed: compute the sum of x^2, y^2 and xz of "posSkin" and build
138 with NumPy the 2x2 matrix
139 inertiaSkin=[[x2,xy], [xy,z2]]
140 Retrieve the eigen vector associated to the maximal eigen value with linalg.eig(inertiaSkin). ::
142 x2=posSkin[:,0]*posSkin[:,0] ; x2=x2.accumulate()[0]
143 y2=posSkin[:,1]*posSkin[:,1] ; y2=y2.accumulate()[0]
144 xy=posSkin[:,0]*posSkin[:,1] ; xy=xy.accumulate()[0]
145 inertiaSkin=matrix([[x2,xy],[xy,y2]])
146 inertiaSkinValues,inertiaSkinVects=linalg.eig(inertiaSkin)
147 pos=max(enumerate(inertiaSkinValues),key=lambda x: x[1])[0]
148 vect0=inertiaSkinVects[pos].tolist()[0]
151 Thanks to the previous computation we can see that the agitator had a rotation of
152 1.1183827931 radian (see solution).
153 Compute and compare the torque on the agitator. ::
155 omega=1.1183827931/(ts[-1][2]-ts[0][2])
156 print("At time-step (%d,%d) at %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))
161 :ref:`python_testmedcouplingloaderex1_solution`