doc/tutorial/medcouplingremapper_fr.rst

   1
   2 MEDCouplingRemapper : interpolation de champs
   3 ---------------------------------------------
   4
   5 Nous allons ici effectuer une interpolation entre deux maillages ``srcMesh`` et ``trgMesh``.
   6 Pour mettre l'accent sur certaines petites subtilités de l'interpolation, nous prenons un cas particulier où ``srcMesh`` est
   7 un maillage raffiné de ``trgMesh`` (avec certaines cellules découpées plus finement).
   8
   9 Début de l'implémentation
  10 ~~~~~~~~~~~~~~~~~~~~~~~~~
  11
  12 Pour commencer l'exercice importer le module Python ``MEDCoupling`` et la
  13 classe ``MEDCouplingRemapper`` du module ``MEDCouplingRemapper``. ::
  14
  15         import MEDCoupling as mc
  16         from MEDCouplingRemapper import MEDCouplingRemapper
  17
  18
  19 Créer le maillage target
  20 ~~~~~~~~~~~~~~~~~~~~~~~~
  21
  22 Construire le maillage non structuré ``trgMesh`` issu d'un maillage cartésien 2D 10x10 commençant au
  23 point ``[0.,0.]`` et ayant un pas de 1. selon X et selon Y : ::
  24
  25         arr = mc.DataArrayDouble(11)
  26         arr.iota(0)
  27         trgMesh = mc.MEDCouplingCMesh()
  28         trgMesh.setCoords(arr,arr)
  29         trgMesh = trgMesh.buildUnstructured()
  30
  31 Créer le maillage source
  32 ~~~~~~~~~~~~~~~~~~~~~~~~
  33
  34 Créer un maillage ``srcMesh`` issu d'un maillage cartésien 2D 20x20 cellules commençant
  35 aussi au point ``[0.,0.]`` et ayant un pas de 0.5 selon X et selon Y : ::
  36
  37         arr = mc.DataArrayDouble(21)
  38         arr.iota(0)
  39         arr *= 0.5
  40         srcMesh = mc.MEDCouplingCMesh()
  41         srcMesh.setCoords(arr,arr)
  42         srcMesh = srcMesh.buildUnstructured()
  43
  44 Afin de rendre l'exercise plus intéressant, triangulariser à l'aide de ``MEDCouplingUMesh.simplexize()``
  45 les 20 premières cellules de ``srcMesh``
  46 (les simplexes 2D sont les triangles). Mettre le résultat dans ``srcMesh``. ::
  47
  48         tmp = srcMesh[:20]    # Extract a sub-part of srcMesh
  49         tmp.simplexize(0)
  50         srcMesh = mc.MEDCouplingUMesh.MergeUMeshes([tmp,srcMesh[20:]])
  51
  52 Interpoler avec MEDCouplingRemapper
  53 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  54
  55 Nous rappelons que pour projeter un champ d'un maillage vers un autre, il faut d'abord préparer la matrice d'interpolation
  56 qui contient les ratios de projection.
  57
  58 Calculer la première partie de la matrice d'interpolation de ``srcMesh`` (discrétisée aux cellules - *P0*)
  59 vers ``trgMesh`` (discrétisée aux cellules également).
  60 Pour ce faire, invoquer ``MEDCouplingRemapper.prepare()`` sur une instance (``remap``) de ``MEDCouplingRemapper``. ::
  61
  62         remap = MEDCouplingRemapper()
  63         remap.prepare(srcMesh,trgMesh,"P0P0")
  64
  65 Vérifier que la matrice calculée par la méthode est correcte dans notre cas trivial.
  66 Pour ce faire, récupérer dans ``myMatrix`` la matrice interne retournée par ``MEDCouplingRemapper.getCrudeMatrix()``.
  67 Celle-ci donne pour chaque cellule de ``trgMesh`` les identifiants de cellules de ``srcMesh`` avec
  68 lesquelles elle s'intersecte, et l'aire d'intersection correspondante.
  69
  70 Vérifier notamment que pour chaque cellule de ``trgMesh`` la somme des aires fait toujours 1. ::
  71
  72         myMatrix = remap.getCrudeMatrix()
  73         print myMatrix
  74         sumByRows = mc.DataArrayDouble(len(myMatrix))
  75         for i,wIt in enumerate(sumByRows):
  76           su = 0.
  77           for it in myMatrix[i]:
  78             su += myMatrix[i][it]
  79           wIt[0] = su
  80         print "Is interpolation well prepared?", sumByRows.isUniform(1.,1e-12)
  81
  82 .. note:: Les triangles dans ``srcMesh`` ont été rajoutés pour casser la monotonie de la matrice ``myMatrix``.
  83
  84 .. note:: Comme on le voit, la préparation ne nécessite que les deux maillages, et rien d'autre.
  85
  86 Construire un champ aux cellules "srcField" construit à partir de la formule analytique
  87 suivante ``7-sqrt((x-5.)*(x-5.)+(y-5.)*(y-5.))`` : ::
  88
  89         srcField = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
  90         srcField.setMesh(srcMesh)
  91         srcField.fillFromAnalytic(1,"7-sqrt((x-5.)*(x-5.)+(y-5.)*(y-5.))")
  92         srcField.getArray().setInfoOnComponent(0, "powercell [W]")
  93
  94 Voici à quoi ressemble ce champ :
  95
  96 .. image:: images/Remapper1.png
  97
  98 Appliquer l'interpolation avec ``MEDCouplingRemapper.transferField()`` : ::
  99
 100         remap.transferField(srcField, 1e300)
 101
 102 .. note:: 1e300 est une valeur par défaut. Cette valeur sera systématiquement assignée à toute cellule
 103         de ``trgField`` n'interceptant aucune cellule de ``srcMesh``. En général, les utilisateurs mettent une
 104         valeur énorme pour repérer ce qui est souvent un bug. Mais d'autres utilisateurs, dans la perspective
 105         d'une interpolation parallèle par exemple, mettent 0.
 106
 107 .. note:: Une exception est envoyée car ``srcField`` n'a pas de *nature* définie.
 108         Nous allons voir dans la suite l'impact de cet attribut sur le résultat final.
 109
 110 Mettre la nature de ``srcField`` à ``IntensiveMaximum``. Cela signifie que le champ doit être interprété commé étant
 111 intensif (une température par exemple). ::
 112
 113         srcField.setNature(mc.IntensiveMaximum)
 114         trgFieldCV = remap.transferField(srcField,1e300)
 115
 116 Vérifier qu'avec la nature ``IntensiveMaximum``, l'intégrale du champ est conservée. Par contre,
 117 la somme sur les cellules (accumulation) n'est **pas** conservée ! ::
 118
 119         integSource = srcField.integral(True)[0]
 120         integTarget =  trgFieldCV.integral(True)[0]
 121         print "IntensiveMaximum -- integrals: %lf == %lf" % (integSource, integTarget)
 122
 123         accSource = srcField.getArray().accumulate()[0]
 124         accTarget = trgFieldCV.getArray().accumulate()[0]
 125         print "IntensiveMaximum -- sums: %lf != %lf" % (accSource, accTarget)
 126
 127
 128 Maintenant mettre la nature de ``srcField`` à ``ExtensiveConservation``. Le champ doit être interprété commé étant
 129 extensif (par exemple une puissance ou un volume). ::
 130
 131         srcField.setNature(mc.ExtensiveConservation)
 132         trgFieldI = remap.transferField(srcField,1e300)
 133
 134 Vérifier qu'avec la nature ``ExtensiveConservation``, l'intégrale du champ n'est **pas** conservée.
 135 Par contre, la somme sur les cellules est conservée. ::
 136
 137         integSource = srcField.integral(True)[0]
 138         integTarget =  trgFieldI.integral(True)[0]
 139         print "ExtensiveConservation -- integrals: %lf != %lf" % (integSource, integTarget)
 140
 141         accSource = srcField.getArray().accumulate()[0]
 142         accTarget = trgFieldI.getArray().accumulate()[0]
 143         print "ExtensiveConservation -- sums: %lf == %lf" % (accSource, accTarget)
 144
 145 Visualiser les champs avec ParaViS, ou en les écrivant dans un fichier.
 146
 147 Solution
 148 ~~~~~~~~
 149
 150 :ref:`python_testMEDCouplingremapper1_solution`