]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/medcoupling_fielddouble1_fr.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / doc / tutorial / medcoupling_fielddouble1_fr.rst
1
2 Manipuler des champs de double
3 ------------------------------
4
5 Les champs dans MEDCoupling ont comme support un unique maillage, de dimension fixée, et bien défini. 
6 Cela semble trivial mais c'est en fait une différence majeure avec la notion de champ dans MED fichier, qui elle est beaucoup
7 plus permissive.
8
9 Les champs sont utiles pour :
10
11 * stocker des valeurs d'une grandeur physique relative au problème traité, mais aussi
12 * des services de haut niveau où l'interaction avec les maillages
13   est requise comme par exemple ``getValueOn()``, ``getValueOnMulti()``, ``integral()``, ``getMeasureField`` 
14   ``normL1()``, ``normL2()``, ``fillFromAnalytic()``, ... qui calculent toutes des valeurs en lien avec le maillage
15   (par exemple le *volume* des cellules)
16 * expliciter précisément les informations échangées entre les différents codes
17   lors de couplage.
18
19 Pour information, l'implémentation de ``MEDCouplingFieldDouble`` est relativement petite car cette classe 
20 délègue la très large majorité de ses traitements à des instances de classes aggrégées 
21 comme ``MEDCouplingMesh``, ``DataArrayDouble``, et ``MEDCouplingSpatialDiscretization``.
22 La classe ``MEDCouplingFieldDouble`` permet d'assurer la cohérence entre tous ces éléments.
23
24
25 Il est souvent  possible et même parfois recommandé de manipuler les tableaux (un ``DataArrayDouble``) 
26 et/ou le maillage d'une instance de ``MEDCouplingFieldDouble`` directement.
27
28 Objectifs
29 ~~~~~~~~~
30
31 Cet exercice met l'accent sur la relation entre le maillage et les valeurs d'un champ.
32
33 * Créer un champ
34 * Agréger des champs
35 * Construire une sous-partie d'un champ
36 * Renuméroter les entités d'un champ
37 * Comparer 2 champs venant de 2 sources différentes
38 * Evaluation d'un champ sur un ensemble de points
39 * Exploser un champ
40
41 Début de l'implementation
42 ~~~~~~~~~~~~~~~~~~~~~~~~~
43
44 Importer le module Python ``MEDCoupling``. ::
45
46         import MEDCoupling as mc
47
48 Créer un ``MEDCouplingUMesh`` à partir d'un maillage 3D cartésien. Chaque direction contiendra 10 cells 
49 et 11 nodes. Le ``MEDCouplingUMesh`` résultant contiendra ainsi 1000 cells. ::
50
51         xarr = mc.DataArrayDouble.New(11,1)
52         xarr.iota(0.)                       # Generate s, s+1, s+2, ... with a given start value s
53         cmesh = mc.MEDCouplingCMesh.New()
54         cmesh.setCoords(xarr,xarr,xarr)
55         mesh = cmesh.buildUnstructured()
56
57 .. note:: La méthode ``MEDCouplingMesh.buildUnstructured()`` est très utile pour construire rapidement un maillage
58         non structuré afin de tester quelque chose.
59
60 Afin de mettre en évidence le problème des types géométriques multiples, convertir en polyhèdres 
61 les cellules d'identifiant pair ::
62
63         mesh.convertToPolyTypes(mc.DataArrayInt.Range(0,mesh.getNumberOfCells(),2))
64
65
66 Création d'un champ
67 ~~~~~~~~~~~~~~~~~~~
68
69 Créer un champ scalaire (une seule composante) aux cellules (P0) appelé "MyField" en appliquant la fonction analytique
70 suivante ``(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)``, où ``(x, y, z)`` représente implicitement les coordonnées du barycentre
71 d'une cellule.
72 Pour cela, deux possiblités :
73
74 * Directement en appelant ``fillFromAnalytic()`` sur un maillage ::
75
76         f = mesh.fillFromAnalytic(mc.ON_CELLS,1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")  # 1 means that the field should have one component
77         f.setName("MyField")
78
79 * Ou en créant au préalable un champ non initialisé, et en appliquant ``fillFromAnalytic()`` sur cette 
80   instance de ``MEDCouplingFieldDouble`` ::
81
82         f2 = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
83         f2.setMesh(mesh)
84         f2.setName("MyField2")
85         f2.fillFromAnalytic(1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)")    # 1 means that the field should have one component
86
87 Comparer les deux champs : comparer ``f`` et ``f2`` avec une précision de 1e-12 sur les coordonnées et
88 de 1e-13 sur les valeurs. ::
89
90         print "Are f and f2 equal?", f.isEqualWithoutConsideringStr(f2,1e-12,1e-13)
91
92
93 .. note:: Le ``WithoutConsideringStr`` dans le nom de la méthode précédente indique que les noms des champs ne seront 
94         pas comparés. On retrouve ce suffixe dans d'autres méthodes MEDCoupling.
95  
96
97 Construire une sous-partie d'un champ
98 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
99
100 Récupérer dans une variable ``ids1`` la liste des identifiants de cellules pour lesquelles la valeur du champ est dans le
101 range [0.0,5.0]. Utiliser pour cela la méthode ``DataArrayDouble.findIdsInRange()``. Avec ce résultat, construire la
102 sous-partie ``fPart1`` du champ ``f``. ::
103
104         da1 = f.getArray()              # a DataArrayDouble, which is a direct reference (not a copy) of the field's values 
105         ids1 = da1.findIdsInRange(0., 5.)
106         fPart1 = f.buildSubPart(ids1)
107         fPart1.writeVTK("ExoField_fPart1.vtu")
108
109 .. image:: images/FieldDouble1_1.png
110
111 Sélectionner la partie ``fPart2`` du champ ``f`` dont toutes les valeurs de tuples
112 sont dans ``[50.,+infinity)``. ::
113
114         ids2 = f.getArray().findIdsInRange(50., 1.e300)
115         fPart2 = f.buildSubPart(ids2)
116
117 Ce genre de technique permet d'extraire facilement les parties d'un champ relatives à un groupe de mailles par exemple.
118
119 Renuméroter les entités d'un champ
120 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
121
122 La partie ``fPart1`` générée est valide d'un point de vue de MEDCoupling. Mais elle 
123 n'est pas valide d'un point de vue de MED *fichier*. 
124 Une renumérotation s'impose dans l'hypothèse de stocker ce champs dans un fichier MED afin d'ordonner les cellules
125 par type géométrique.
126
127 L'idée est d'utiliser les deux méthodes ``MEDCouplingUMesh.sortCellsInMEDFileFrmt()`` et
128 ``DataArrayDouble.renumberInPlace()`` pour renuméroter manuellement une *copie* de ``fPart1``: ::
129
130         fPart1Cpy = fPart1.deepCopy()
131         o2n = fPart1Cpy.getMesh().sortCellsInMEDFileFrmt()
132         fPart1Cpy.getArray().renumberInPlace(o2n)
133
134 ``fPart1Cpy`` est désormais normalisé pour être stocké dans un fichier MED (ce que nous verrons plus loin)
135
136 Vérifier que ``fPart1Cpy`` et ``fPart1`` sont les mêmes à une permutation près (``MEDCouplingFieldDouble.substractInPlaceDM()``) ::
137
138         fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
139         fPart1Cpy.getArray().abs()
140         print "Equal field ? %s" % (fPart1Cpy.getArray().accumulate()[0]<1e-12)
141
142 .. note:: La renumérotation effectuée ici représente en fait d'un cas très particulier
143         d'interpolation. Effectivement l'hypothèse est faite que les supports
144         de ``fPart1`` et ``fPart1Cpy`` sont égaux à une permutation de cellule
145         et/ou noeuds.  
146
147 Agréger des champs
148 ~~~~~~~~~~~~~~~~~~
149
150 Agréger ``fPart1`` et ``fPart2`` (utiliser ``MEDCouplingFieldDouble.MergeFields()``). Et mettre le résultat de l'agrégation
151 dans ``fPart12``. ::
152
153         fPart12 = mc.MEDCouplingFieldDouble.MergeFields([fPart1,fPart2])
154         fPart12.writeVTK("ExoField_fPart12.vtu")
155
156 .. note:: La méthode ``MEDCouplingFieldDouble.MergeFields()`` devrait vraiment se 
157         nommer ``MEDCouplingFieldDouble.AggregateFields()`` ...
158
159 .. image:: images/FieldDouble1_2.png
160
161 Evaluation d'un champ en des points donnés de l'espace
162 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
163
164 Evaluer la valeur du champ ``fPart12`` calculé précédemment sur les barycentres des cellules de son
165 maillage (variable ``bary``) et mettre le résultat dans ``arr1``.
166 Utiliser pour cela les méthodes ``MEDCouplingFieldDouble.getValueOnMulti()`` et ``MEDCouplingMesh.computeCellCenterOfMass()``.  
167
168 De manière similaire, évaluer ensuite directement le champ ``f`` en utilisant la même liste de points
169 que précédemment (``bary``) et mettre le résultat dans ``arr2``.
170
171 Vérifier ensuite que ``arr1`` et ``arr2`` sont bien égaux: ::
172
173         bary = fPart12.getMesh().computeCellCenterOfMass()
174         arr1 = fPart12.getValueOnMulti(bary)
175         arr2 = f.getValueOnMulti(bary)
176         delta = arr1-arr2
177         delta.abs()
178         print "Is field evaluation matching?", (delta.accumulate()[0]<1e-12)
179
180 .. note:: Dans ce contexte, et pour un champ aux cellules (P0) par exemple, "évaluer" en un point signifie retourner la valeur 
181         de la cellule contenant le point donné.
182         Pour les champs aux noeuds (P1), les cellules doivent être de types simples (triangles, tétraèdres) et une interpolation
183         linéaire est alors utilisée.
184
185 .. note:: Cette technique peut être utilisée pour juger rapidement de la qualité d'une interpolation.
186
187 Opérations sur les champs
188 ~~~~~~~~~~~~~~~~~~~~~~~~~
189
190 Calculer l'intégrale du champ ``fPart12`` sur le maillage, et la retrouver d'une autre manière en utilisant
191 la méthode ``DataArrayDouble.accumulate()`` sur le tableau de valeurs de ce champ. 
192 On rappelle que, vu le maillage simplifié en jeu, les cellules ont toutes un volume unité. :: 
193
194         integ1 = fPart12.integral(0,True)
195         integ1_bis = fPart12.getArray().accumulate()[0]
196         print "First integral matching ?", ( abs(integ1 - integ1_bis) < 1e-8 )
197
198 Ensuite appliquer une homotétie de facteur 1.2 centrée en [0.,0.,0.] sur le support de ``fPart12`` (c'est-à-dire son maillage).
199 Quelle est alors la nouvelle valeur de l'intégrale ? ::
200
201         fPart12.getMesh().scale([0.,0.,0.], 1.2)
202         integ2 = fPart12.integral(0,True)
203         print "Second integral matching ?", ( abs(integ2-integ1_bis*1.2*1.2*1.2) < 1e-8 )
204
205 Exploser un champ - Vecteurs de déplacement
206 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
207
208 Nous allons maintenant créer un nouveau maillage représentant l'*éclaté* du maillage initial.
209
210 Partant du maillage ``mesh`` créer un champ vectoriel aux cellules ``fVec`` ayant 3 composantes représentant 
211 le vecteur déplacement entre le point [5.,5.,5.] et le barycentre de chaque cellule du maillage.
212 Utiliser la méthode ``MEDCouplingMesh.fillFromAnalytic()`` : ::
213
214         fVec = mesh.fillFromAnalytic(mc.ON_CELLS,3,"(x-5.)*IVec+(y-5.)*JVec+(z-5.)*KVec")
215
216 .. note:: Les identifiants spéciaux ``IVec``, ``JVec`` et ``KVec`` représentent les vecteurs unitaires du repère. 
217
218 Créer ensuite une réduction de ``fVec`` (nommée ``fVecPart1``) sur les cellules ``ids1`` précédemment obtenues : ::
219
220         fVecPart1 = fVec.buildSubPart(ids1)
221         fVecPart1.setName("fVecPart1")
222
223 Construire le champ scalaire ``fPart1Exploded`` ayant les mêmes valeurs que ``fPart1`` mais reposant sur un maillage *eclaté*
224 par rapport à celui de ``fPart1.getMesh()``. Pour exploser ``fPart1.getMesh()`` utiliser le champ de déplacement vectoriel
225 ``fVecPart1`` afin d'appliquer à chaque cellule la translation associée. ::
226
227         cells = fPart1.getMesh().getNumberOfCells() * [None]
228         for icell,vec in enumerate(fVecPart1.getArray()):
229           m = fPart1.getMesh()[[icell]]
230           m.zipCoords()      # Not mandatory but saves memory
231           m.translate(vec)
232           cells[icell] = m
233           pass
234         meshFVecPart1Exploded = mc.MEDCouplingUMesh.MergeUMeshes(cells)
235         fPart1.setMesh(meshFVecPart1Exploded)
236         fPart1.writeVTK("ExoField_fPart1_explo.vtu")
237
238 Et voilà ce que vous devriez obtenir:
239
240 .. image:: images/FieldDouble1_1_exploded.png
241         :scale: 120
242         
243 Solution
244 ~~~~~~~~
245
246 :ref:`python_testMEDCouplingfielddouble1_solution`