doc/tutorial/medcouplingloaderex1_en.rst

   1
   2
   3 Full example 1 - Agitator
   4 -------------------------
   5
   6 The MED file :download:`agitateur.med <data/agitateur.med>` is used and has the following content:ant :
   7
   8 .. image:: images/agitateur.jpg
   9
  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
  12 to the other.
  13 The purpose of the exercise is to compute the torque applied on this piece.agitateur.
  14
  15 Objective
  16 ~~~~~~~~~
  17
  18 The aim of this exercise is to give a full example of non-trivial post-treatment
  19 from a MED file.
  20
  21 Implementation start
  22 ~~~~~~~~~~~~~~~~~~~~
  23
  24 Import the whole Python module MEDLoader (which includes MEDCoupling).
  25 Also import NumPy and acos() from the math module. ::
  26
  27         from MEDLoader import *
  28         from numpy import *
  29         from math import acos
  30
  31 Mesh and field extraction using advanced API
  32 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  33
  34 Using the advanced API read the whole file "agitateur.med" and display all time-steps of
  35 the first field. ::
  36
  37         data=MEDFileData("agitateur.med")
  38         ts=data.getFields()[0].getTimeSteps()
  39         print ts
  40
  41 Get the agitator's mesh (in green) at the time-step (2,-1) (see ts).
  42 To this end use the cell field "DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM" and select
  43 only the field part having a value within [0.0, 1.0] (variable "ids"). ::
  44
  45         fMts=data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
  46         f1ts=fMts[(2,-1)]
  47         fMc=f1ts.getFieldAtLevel(ON_CELLS,0)
  48         arr=fMc.getArray()
  49         arr.getMinMaxPerComponent() # just to see the variation range of the field per component
  50         ids=arr.getIdsInRange(0.,1.)
  51         f2Mc=fMc[ids]
  52
  53 Using the field "PRESSION_ELEM_DOM" find the 3D pression field applied on the agitator.
  54 Store the result in pressOnAgitateur. ::
  55
  56         pressMts=data.getFields()["PRESSION_ELEM_DOM"]
  57         press1ts=pressMts[(2,-1)]
  58         pressMc=press1ts.getFieldAtLevel(ON_CELLS,0)
  59         pressOnAgitateurMc=pressMc[ids]
  60
  61 Delete unused nodes in pressOnAgitateurMc.getMesh(). ::
  62
  63         pressOnAgitateurMc.getMesh().zipCoords()
  64
  65 Create a 3D surface field from the 3D cell field
  66 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  67
  68 Deduce the 3D field on the skin of the agitator.
  69 To achieve this use the constituting mesh MEDCouplingUMesh.buildDescendingConnectivity().
  70 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  71
  72         agitateurMesh3DMc=pressOnAgitateurMc.getMesh()
  73         m3DSurf,desc,descI,revDesc,revDescI=agitateurMesh3DMc.buildDescendingConnectivity()
  74         nbOf3DCellSharing=revDescI.deltaShiftIndex()
  75         ids2=nbOf3DCellSharing.getIdsEqual(1)
  76         agitateurSkinMc=m3DSurf[ids2]
  77         OffsetsOfTupleIdsInField=revDescI[ids2]
  78         tupleIdsInField=revDesc[OffsetsOfTupleIdsInField]
  79         pressOnSkinAgitateurMc=pressOnAgitateurMc[tupleIdsInField]
  80         pressOnSkinAgitateurMc.setMesh(agitateurSkinMc)
  81
  82 Manipulate fields
  83 ~~~~~~~~~~~~~~~~~
  84
  85 Compute the force vector field on the agitator's skin by multiplying at each cell
  86 the pressure by the surface, and then the normal vector.
  87 Pression is expressed in bar, convert it first to Pa. ::
  88
  89         pressSkin=pressOnSkinAgitateurMc.getArray()
  90         pressSkin*=1e5
  91         areaSkin=agitateurSkinMc.getMeasureField(True).getArray()
  92         forceSkin=pressSkin*areaSkin
  93         normalSkin=agitateurSkinMc.buildOrthogonalField().getArray()
  94         forceVectSkin=forceSkin*normalSkin
  95
  96 First computation of the torque at the center of mass of the agitator:
  97
  98 Let's compute first the position of the center of mass.
  99 Compute the polyhedron representing the 3D mesh hull of the agitator "agitateurMesh3DMc"
 100 (use MEDCouplingUMesh.buildSpreadZonesWithPoly()). ::
 101
 102         singlePolyhedron=agitateurMesh3DMc.buildSpreadZonesWithPoly()
 103         singlePolyhedron.orientCorrectlyPolyhedrons()
 104         centerOfMass=singlePolyhedron.getBarycenterAndOwner()
 105
 106 .. note:: The call to MEDCouplingUMesh.orientCorrectlyPolyhedrons() is not mandatory
 107         but is recommended: if the polyhedron happens to be mis-oriented, its center of mass will
 108         be incorrect!
 109
 110 Compute for each skin cell  the torque with respect to the center of mass "centerOfMass".
 111 To this end compute "posSkin", a DataArrayDouble giving for each skin cell the vector
 112 centerOfMass -> G, where G represents the center of mass of the current cell. ::
 113
 114         barySkin=agitateurSkinMc.getBarycenterAndOwner()
 115         posSkin=barySkin-centerOfMass
 116
 117 Compute the cross product for each cell of "posSkin" using "forceVectSkin"
 118 (method DataArrayDouble.CrossProduct()). ::
 119
 120         torquePerCellOnSkin=DataArrayDouble.CrossProduct(posSkin,forceVectSkin)
 121
 122 Sum "torqueOnSkin" using DataArrayDouble.accumulate(). ::
 123
 124        zeTorque=torquePerCellOnSkin.accumulate()
 125        print "couple = %r N.m"%(zeTorque[2])
 126
 127 Check the previously computed torque by dividing the power by the angular speed.
 128 Compute the power per skin cell and sum it. ::
 129
 130        speedMts=data.getFields()["VITESSE_ELEM_DOM"]
 131        speed1ts=speedMts[(2,-1)]
 132        speedMc=speed1ts.getFieldAtLevel(ON_CELLS,0)
 133        speedOnSkin=speedMc.getArray()[tupleIdsInField]
 134        powerSkin=DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
 135        power=powerSkin.accumulate()[0]
 136        print "power = %r W"%(power)
 137
 138 Compute the angular speed: compute the sum of x^2, y^2 and xz of "posSkin" and build
 139 with NumPy the 2x2 matrix
 140 inertiaSkin=[[x2,xy], [xy,z2]]
 141 Retrieve the eigen vector associated to the maximal eigen value with linalg.eig(inertiaSkin). ::
 142
 143        x2=posSkin[:,0]*posSkin[:,0] ; x2=x2.accumulate()[0]
 144        y2=posSkin[:,1]*posSkin[:,1] ; y2=y2.accumulate()[0]
 145        xy=posSkin[:,0]*posSkin[:,1] ; xy=xy.accumulate()[0]
 146        inertiaSkin=matrix([[x2,xy],[xy,y2]])
 147        inertiaSkinValues,inertiaSkinVects=linalg.eig(inertiaSkin)
 148        pos=max(enumerate(inertiaSkinValues),key=lambda x: x[1])[0]
 149        vect0=inertiaSkinVects[pos].tolist()[0]
 150        print vect0
 151
 152 Thanks to the previous computation we can see that the agitator had a rotation of
 153 1.1183827931 radian (see solution).
 154 Compute and compare the torque on the agitator. ::
 155
 156        omega=1.1183827931/(ts[-1][2]-ts[0][2])
 157        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)
 158
 159 Solution
 160 ~~~~~~~~
 161
 162 :ref:`python_testmedcouplingloaderex1_solution`