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