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