]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/medcoupling_fielddouble1_en.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / doc / tutorial / medcoupling_fielddouble1_en.rst
1
2 Playing with fields
3 -------------------
4
5 In MEDCoupling fields are directly supported by a mesh. That is a major difference
6 with the field notion in MED file.
7
8 Fields instances are well suited:
9   * for high level services where an interaction with the mesh is requested as in getValueOn(), getValueOnMulti(), integral(), normL1(), normL2(), fillFromAnalytic(), changeUnderlyingMesh(), ...
10   * to precisely transmit information between components when coupling codes.
11
12 For information, the MEDCouplingFieldDouble implementation is tiny since 
13 it delegates most of its work to MEDCouplingMesh, DataArrayDouble, MEDCouplingSpatialDiscretization classes.
14 MEDCouplingFieldDouble insure the coherency in this model.
15 It is often possible and even advised to manipulate the array and mesh of a MEDCouplingFieldDouble directly.
16
17 Objective
18 ~~~~~~~~~
19
20 This exercise focuses on the relationship between a mesh and the values of a field.
21
22 * Create a field
23 * Aggregate fields
24 * Build parts of a field
25 * Renumbering a field
26 * Compare two fields coming from 2 sources
27 * Evaluation of a field on a point set
28 * Explode a field 
29
30 Implementation start
31 ~~~~~~~~~~~~~~~~~~~~
32
33 Import the MEDCoupling Python module. ::
34
35         from MEDCoupling import *
36
37 We are going to create a MEDCouplingUMesh from a 3D cartesian mesh. Each direction will contain 10 cells and 11 nodes. The generated MEDCouplingUMesh
38 will contain 1000 cells. ::
39
40         xarr=DataArrayDouble.New(11,1)
41         xarr.iota(0.)
42         cmesh=MEDCouplingCMesh.New()
43         cmesh.setCoords(xarr,xarr,xarr)
44         mesh=cmesh.buildUnstructured()
45         
46 In order to put the focus on mixed types, cells with an even id will be converted to polyhedrons ::
47
48         mesh.convertToPolyTypes(DataArrayInt.Range(0,mesh.getNumberOfCells(),2))
49         
50 Creation of a field
51 ~~~~~~~~~~~~~~~~~~~
52
53 Create a field called "MyField" by application of an analytic function "(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)".
54 Two possibilities :
55
56 * Directly by calling fillFromAnalytic on a mesh ::
57
58         f=mesh.fillFromAnalytic(ON_CELLS,1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")
59         f.setName("MyField")
60         
61 * Or by building a new instance ::
62
63         f2=MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME)
64         f2.setMesh(mesh)
65         f2.setName("MyField2")
66         f2.fillFromAnalytic(1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")
67
68 Compare the two fields:
69 Compare f and f2 with a precision of 1e-12 on coordinates and 1e-12 on values. ::
70
71         print "f and f2 are equal: %s"%(f.isEqualWithoutConsideringStr(f2,1e-12,1e-12))
72
73 Builing of a subpart of a field
74 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
75         
76 Store in ids1 the list of tuple ids whose value is within [0.0,5.0] (DataArrayDouble.findIdsInRange)    . From ids1 build the sub-part fPart1 of the field "f". ::
77
78         ids1=f.getArray().findIdsInRange(0.,5.)
79         fPart1=f.buildSubPart(ids1)
80         
81 .. image:: images/FieldDouble1_1.png
82
83 Select the part "fPart2" of the field "f" whose values are in [50.,infinity). ::
84
85         ids2=f.getArray().findIdsInRange(50.,1.e300)
86         fPart2=f.buildSubPart(ids2)
87
88 Renumbering cells
89 ~~~~~~~~~~~~~~~~~
90
91 The generated file "fPart1" is valid for MEDCoupling, but its cells are not sorted by geometric type: it is not valid from a MED file point of view. By using MEDCouplingUMesh.sortCellsInMEDFileFrmt and DataArrayDouble.renumberInPlace
92 renumber manually fPart1 starting from a deep copy of fPart1. ::
93
94         fPart1Cpy=fPart1.deepCopy()
95         o2n=fPart1Cpy.getMesh().sortCellsInMEDFileFrmt()
96         fPart1Cpy.getArray().renumberInPlace(o2n)
97         
98 "fPart1Cpy" is now normalized to be stored in a MED file (we will tackle this later).   
99 Check that fPart1Cpy and fPart1 are the same (discarding any permutation): ::
100
101         fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
102         fPart1Cpy.getArray().abs()
103         print "Fields are the same? %s"%(fPart1Cpy.getArray().accumulate()[0]<1e-12)
104
105 .. note:: This is in fact a very special case of interpolation. Except that here 
106         we assume that the supports of "fPart1" and "fPart1Cpy" are equal, discarding any 
107         cell and/or node permutation.
108
109 Aggregate Fields
110 ~~~~~~~~~~~~~~~~
111
112 Aggregate fields "fPart1" and "fPart2". The result is stored in "fPart12". ::
113
114         fPart12=MEDCouplingFieldDouble.MergeFields([fPart1,fPart2])
115
116 .. image:: images/FieldDouble1_2.png
117
118 .. note:: Apologies for the name MEDCouplingFieldDouble.MergeFields instead of 
119                 AggregateFields.
120
121 Evaluation of a MEDCouplingFieldDouble on given space points
122 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
123
124 Evaluate the values of the computed field "fPart12" on the barycenters of its mesh.
125 Evaluate the field "f" on the same barycenters. The method used is MEDCouplingFieldDouble.getValueOnMulti(). ::
126
127         bary=fPart12.getMesh().computeCellCenterOfMass()
128         arr1=fPart12.getValueOnMulti(bary)
129         arr2=f.getValueOnMulti(bary)
130         delta=arr1-arr2
131         delta.abs()
132         print "Check OK: %s"%(delta.accumulate()[0]<1e-12)
133
134 .. note:: In this context and for example for a field on cells, "evaluate" at a point means returning the value of the cell containing the point.
135 .. note:: This technique can be used to quickly assess the quality of an interpolation.
136
137 Operations on a field
138 ~~~~~~~~~~~~~~~~~~~~~
139
140 Compute the integral of the field "fPart12" and compute it a second time by using
141 DataArrayDouble.accumulate on the underlying DataArrayDouble of this "fPart12" (remember that the cell volumes are all 1.0). 
142 To show the link with the underlying mesh, scale the underlying mesh (fPart12.getMesh()) by 1.2 and centered at [0.,0.,0.].
143 Recompute the integral.
144 ::
145
146         fPart12.integral(0,True)
147         fPart12.getArray().accumulate()
148         fPart12.getMesh().scale([0.,0.,0.],1.2)
149         abs(fPart12.integral(0,True)-fPart12.getArray().accumulate()[0]*1.2*1.2*1.2)<1e-8
150
151 Exploding a field
152 ~~~~~~~~~~~~~~~~~
153
154 Starting from "mesh", create a vector field on cells "fVec" with 3 components representing the displacement between each cell's barycenter and the point [5.,5.,5.]. Use MEDCouplingMesh.fillFromAnalytic(). ::
155
156         fVec=mesh.fillFromAnalytic(ON_CELLS,3,"(x-5.)*IVec+(y-5.)*JVec+(z-5.)*KVec")
157
158 Create the reduction of "fVec" ("fVecPart1") on cell IDs "ids1" (previously obtained). ::
159
160         fVecPart1=fVec.buildSubPart(ids1)
161         fVecPart1.setName("fVecPart1")
162
163 Build the scalar field fPart1Exploded having the same values as "fPart1" but supported by an exploded mesh (in comparison to fPart1.getMesh()).
164 To explode the underlying mesh fPart1.getMesh(), use the vectorial displacement field "fVecPart1" in order to apply to each cell the proper translation. ::
165
166         cells=fPart1.getMesh().getNumberOfCells()*[None]
167         for icell,vec in enumerate(fVecPart1.getArray()):
168           m=fPart1.getMesh()[[icell]]
169           m.zipCoords() # pas absolument nécessaire mais permet d'être économe en mémoire
170           m.translate(vec)
171           cells[icell]=m
172           pass
173         meshFVecPart1Exploded=MEDCouplingUMesh.MergeUMeshes(cells)
174         fPart1.setMesh(meshFVecPart1Exploded)
175
176 .. image:: images/FieldDouble1_1_exploded.png
177         
178 Solution
179 ~~~~~~~~
180
181 :ref:`python_testMEDCouplingfielddouble1_solution`