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 medcoupling.
  25 Also import NumPy. ::
  26
  27     import medcoupling as mc
  28     import numpy as np
  29
  30 Mesh and field extraction using advanced API
  31 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  32
  33 Using the advanced API read the whole file "agitateur.med" and display all time-steps of
  34 the first field. ::
  35
  36         data=mc.MEDFileData("agitateur.med")
  37         ts=data.getFields()[0].getTimeSteps()
  38         print(ts)
  39
  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"). ::
  43
  44         fMts=data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
  45         f1ts=fMts[(2,-1)]
  46         fMc=f1ts.getFieldAtLevel(mc.ON_CELLS,0)
  47         arr=fMc.getArray()
  48         arr.getMinMaxPerComponent() # just to see the variation range of the field per component
  49         ids=arr.findIdsInRange(0.,1.)
  50         f2Mc=fMc[ids]
  51
  52 Using the field "PRESSION_ELEM_DOM" find the 3D pression field applied on the agitator.
  53 Store the result in pressOnAgitateur. ::
  54
  55         pressMts=data.getFields()["PRESSION_ELEM_DOM"]
  56         press1ts=pressMts[(2,-1)]
  57         pressMc=press1ts.getFieldAtLevel(mc.ON_CELLS,0)
  58         pressOnAgitateurMc=pressMc[ids]
  59
  60 Delete unused nodes in pressOnAgitateurMc.getMesh(). ::
  61
  62         pressOnAgitateurMc.getMesh().zipCoords()
  63
  64 Create a 3D surface field from the 3D cell field
  65 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  66
  67 Deduce the 3D field on the skin of the agitator.
  68 To achieve this use the constituting mesh MEDCouplingUMesh.buildDescendingConnectivity().
  69 :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  70
  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)
  80
  81 Manipulate fields
  82 ~~~~~~~~~~~~~~~~~
  83
  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. ::
  87
  88         pressSkin=pressOnSkinAgitateurMc.getArray()
  89         pressSkin*=1e5
  90         areaSkin=agitateurSkinMc.getMeasureField(True).getArray()
  91         forceSkin=pressSkin*areaSkin
  92         normalSkin=agitateurSkinMc.buildOrthogonalField().getArray()
  93         forceVectSkin=forceSkin*normalSkin
  94
  95 First computation of the torque at the center of mass of the agitator:
  96
  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()). ::
 100
 101         singlePolyhedron=agitateurMesh3DMc.buildSpreadZonesWithPoly()
 102         singlePolyhedron.orientCorrectlyPolyhedrons()
 103         centerOfMass=singlePolyhedron.computeCellCenterOfMass()
 104
 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
 107         be incorrect!
 108
 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. ::
 112
 113         barySkin=agitateurSkinMc.computeCellCenterOfMass()
 114         posSkin=barySkin-centerOfMass
 115
 116 Compute the cross product for each cell of "posSkin" using "forceVectSkin"
 117 (method DataArrayDouble.CrossProduct()). ::
 118
 119         torquePerCellOnSkin=mc.DataArrayDouble.CrossProduct(posSkin,forceVectSkin)
 120
 121 Sum "torqueOnSkin" using DataArrayDouble.accumulate(). ::
 122
 123        zeTorque=torquePerCellOnSkin.accumulate()
 124        print("couple = %r N.m"%(zeTorque[2]))
 125
 126 Check the previously computed torque by dividing the power by the angular speed.
 127 Compute the power per skin cell and sum it. ::
 128
 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))
 136
 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). ::
 141
 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]
 149        print(vect0)
 150
 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. ::
 154
 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))
 157
 158 Solution
 159 ~~~~~~~~
 160
 161 :ref:`python_testmedcouplingloaderex1_solution`