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`