]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/medcouplingloaderex2_fr.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / doc / tutorial / medcouplingloaderex2_fr.rst
1
2 MEDCoupling / MEDLoader - Exemple complet 2 - RJH
3 -------------------------------------------------
4
5 Ici nous partons de deux fichiers MED très librement inspirés du réacteur expérimental
6 RJH (Réacteur Jules Horowitz).
7
8 Le premier, :download:`Fixe.med <data/Fixe.med>`, représente la partie statique du réacteur
9 sans dispositif expérimental.  
10
11 .. image:: images/fixm.jpg
12         :scale: 70
13
14 Le deuxième, :download:`Mobile.med <data/Mobile.med>`, représente la partie mobile.
15
16 .. image:: images/mobm.jpg
17         :scale: 70
18
19
20 Objectif
21 ~~~~~~~~
22
23 Le but ici est d'utiliser MEDCoupling pour:
24
25 * intersecter ces deux maillages,
26 * y mettre un champ,
27 * et ainsi localiser les zones de recouvrement
28
29
30 Début de l'implémentation
31 ~~~~~~~~~~~~~~~~~~~~~~~~~
32
33 Pour commencer l'exercice importer tout le module python MEDLoader (qui inclus MEDCoupling). ::
34
35         import MEDLoader as ml
36
37 Lire et réparer le maillage statique "Fixe.med"
38 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
39
40 Avec l'API avancée lire tout le fichier "Fixe.med" et appeler ``fixm``
41 l'objet de type ``MEDCouplingUMesh`` représentant le maillage statique. ::
42
43         fixe = ml.MEDFileMesh.New("Fixe.med")
44         fixm = fixe.getMeshAtLevel(0)
45
46 Pour ce qui suit il faut absolument que deux cellules se touchant partagent les mêmes edges. Pour ce faire, comme on est
47 en connectivité nodale, il faut absolument que les noeuds soient les mêmes. Il s'avère que cela n'est pas le cas ici.
48 Fusionner le noeuds distants de moins de 1e-10 et regarder l'impact sur le nombre de noeuds de ``fixm``. ::
49
50         print "Nb of nodes in the file : %i " % (fixm.getNumberOfNodes())
51         fixm.mergeNodes(1e-10)
52         print "Nb of non duplicated nodes : %i" % (fixm.getNumberOfNodes())
53
54 Même traitement pour ``Mobile.med``, le lire avec l'API avancée de MEDLoader (appeler ``mobm`` l'instance du maillage) 
55 et le réparer en supprimant les noeuds dupliqués. ::
56
57         mobile = ml.MEDFileMesh.New("Mobile.med")
58         mobm = mobile.getMeshAtLevel(0)
59         mobm.mergeNodes(1e-10)
60
61 Le maillage du RJH étant plus général que des ``TRI6`` et des ``QUAD8``, on a besoin
62 de stocker ces cellules avec un type géométrique à ``QPOLYG`` (Quadratic Polygon) qui représente un polygone *quadratique* 
63 (le terme n'est pas très heureux, encore des raisons historiques, ...), c'est-à-dire un polygone avec un nombre arbitraire
64 de côtés, et potentiellement des côtés en forme d'arcs de cercles plutôt que de segments de droites.
65 Ce type géométrique ``NORM_QPOLYG`` est dans MEDCoupling/MEDLoader et aussi dans MED fichier.
66
67 Nous voudrions visualiser ces deux maillages dans PARAVIS/ParaView, mais nous rencontrons ici deux soucis:
68
69 * les polygones non-convexes sont, par défaut, mal représentés par VTK en mode *Surface*.
70   Il faut sélectionner l'option avancée "Triangulate" dans le panneau Display de PARAVIS/ParaView pour avoir un rendu correct.
71 * les arcs de cercles ne sont pas correctement supportés par ParaView. Il faut les *tesséliser*, c'est-à-dire les transformer
72   en plusieurs petits segments de droite. La méthode ``MEDCouplingUMesh.tessellate2D()`` fait ce travail, mais modifie
73   le maillage. Nous faisons donc une copie préalable. Nous passons en paramètre la finesse de découpage (0.1 est suffisant 
74   -- angle en radian). Attention donc à ne pas modifer ni ``fixm`` ni ``mobm`` ! ::
75
76         fixm2 = fixm.deepCopy()        # tessellate2D() modifies the current mesh
77         fixm2.tessellate2D(0.1)
78         fixm2.writeVTK("fixm2.vtu")
79         mobm2 = mobm.deepCopy()
80         mobm2.tessellate2D(0.1)
81         mobm2.writeVTK("mobm2.vtu")
82
83 Faire une petite méthode ``displayVTK()``, faisant le travail qui nous servira souvent après. ::
84
85         def displayVTK(m,fname):
86                 tmp = m.deepCopy()
87                 tmp.tessellate2D(0.1)
88                 tmp.writeVTK(fname)
89                 return
90
91 Faire des réductions et des repérages de zones
92 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
93
94 Le maillage ``mobm`` est en 6 parties distinctes (voir l'image au dessus). On ne veut récupérer que la première partie.
95 Utiliser ``MEDCouplingUMesh.partitionBySpreadZone()`` pour partitionner en zones ``mobm`` et ne prendre que la première zone.
96 Appeler cette nouvelle instance ``zone1Mobm`` et lui retirer tous les noeuds orphelins (``MEDCouplingUMesh.zipCoords()``) 
97 Enfin l'afficher : ::
98
99         zonesInMobm = mobm.partitionBySpreadZone()
100         print "Nb of zones in mobm : %i" % (len(zonesInMobm))
101         zone1Mobm = mobm[zonesInMobm[0]]
102         zone1Mobm.zipCoords()
103         displayVTK(zone1Mobm, "zone1Mobm.vtu")
104
105 .. image:: images/zone1Mobm.jpg
106         :scale: 70
107
108 Nous allons désormais travailler autour de ``zone1Mobm``. Nous allons réduire la zone de travail de ``fixm`` autour de ``zone1Mobm``.
109 Pour ce faire, réduire ``fixm`` en ne prenant que les cellules dans la boîte englobante 
110 de ``zone1Mobm`` (``MEDCouplingUMesh.getBoundingBox()`` et ``MEDCouplingUMesh.getCellsInBoundingBox()``).
111 Appeler ce nouvel objet ``partFixm``, lui retirer ses noeuds orphelins et l'afficher. ::
112
113         ids2 = fixm.getCellsInBoundingBox(zone1Mobm.getBoundingBox(),1e-10)
114         partFixm = fixm[ids2]
115         partFixm.zipCoords()
116         displayVTK(partFixm,"partFixm.vtu")
117
118 .. image:: images/partFixmAndzone1Mobm.jpg
119
120 Intersecter géométriquement deux maillages
121 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
122
123 C'est le coeur de l'exercice. Nous allons intersecter géométriquement ``partFixm`` et ``zone1Mobm``. Cela revient à 
124 partitionner à minima ``partFixm`` en cellules appartenant
125 soit complètement à ``partFixm`` soit à ``partFixm`` et ``zone1Mobm``. Invoquer la méthode statique 
126 ``MEDCouplingUMesh.Intersect2DMeshes()``, avec ``partFixm`` et ``zone1Mobm`` et mettre une précision
127 de 1e-10 (seuil de détection de fusion). Cette méthode retourne 3 paramètres (voir API dans la doc) que l'on appellera 
128 ici ``partFixMob``, ``iPart`` et ``iMob`` dans cet ordre.
129
130 Sur ``partFixMob`` merger les noeuds à 1e-10 près. ::
131
132         partFixMob, iPart, iMob = ml.MEDCouplingUMesh.Intersect2DMeshes(partFixm,zone1Mobm,1e-10)
133         partFixMob.mergeNodes(1e-10)
134
135 Récupérer et afficher la partie de ``partFixm`` qui n'est pas dans ``zone1Mobm``. Appeler ce maillage ``partFixmWithoutZone1Mobm``. ::
136
137         ids3 = iMob.findIdsEqual(-1)
138         partFixmWithoutZone1Mobm = partFixMob[ids3]
139         displayVTK(partFixmWithoutZone1Mobm,"partFixmWithoutZone1Mobm.vtu")
140
141 .. image:: images/partFixmWithoutZone1Mobm.jpg
142         :scale: 70
143
144 Maintenant, on va vérifier la qualité du résultat retourné par ``MEDCouplingUMesh.Intersect2DMeshes()``. 
145 Pour ce faire on va passer 3 tests:
146
147  * **Check #0** la somme des aires des cellules de ``partFixm`` et égale à celle de ``partFixMob``
148  * **Check #1** la somme des aires des cellules de ``zone1Mobm`` et égale à la somme des cells de ``partFixMob`` 
149    dont l'id dans ``iMob`` est different de -1
150  * **Check #2** pour chaque cellule de ``partFixm``, son aire est égale à la somme des aires des cellules de ``partFixMob``
151
152 L'aire est une valeur algébrique. Donc attention cette verification ne peut se faire que si les cellules 
153 sont toutes bien orientées ou à minima toutes orientées de la même manière.
154 Pour ce faire, regardons les aires des 38 cellules de ``partFixm`` (nom de variable : ``areaPartFixm``). ::
155
156         areaPartFixm = partFixm.getMeasureField(ml.ON_CELLS).getArray()
157         print areaPartFixm.getValues()
158
159 On voit que toutes les valeurs sont négatives. *Bilan*: ce fichier MED ne respecte pas la convention MED fichier !
160 ``partFixm`` étant mal orienté, et ``MEDCouplingUMesh.Intersect2DMeshes()`` conservant l'orientation, 
161 ``partFixMob`` est lui aussi mal orienté.
162 Bref, on va faire les comparaisons sur des tableaux de valeurs absolues. Vérifier alors **Check #0**. ::
163
164         areaPartFixm = partFixm.getMeasureField(ml.ON_CELLS).getArray()
165         areaPartFixm.abs()
166         areaPartFixMob = partFixMob.getMeasureField(ml.ON_CELLS).getArray()
167         areaPartFixMob.abs()
168         val1=areaPartFixm.accumulate()[0]
169         val2=areaPartFixMob.accumulate()[0]
170         print "Check #0 %lf == %lf with precision 1e-8? %s" % (val1,val2,str(abs(val1-val2)<1e-8))
171
172 On peut passer au **Check #1**. L'esprit est le même que le **Check #0**. ::
173
174         areaZone1Mobm = zone1Mobm.getMeasureField(ml.ON_CELLS).getArray()
175         areaZone1Mobm.abs()
176         val3 = areaZone1Mobm.accumulate()[0]
177         ids4 = iMob.findIdsNotEqual(-1)
178         areaPartFixMob2 = areaPartFixMob[ids4]
179         val4 = areaPartFixMob2.accumulate()[0]
180         print "Check #1 %lf == %lf with precision 1e-8 ? %s" % (val3,val4,str(abs(val3-val4)<1e-8))
181
182 Puis le **Check #2**. ::
183
184         isCheck2OK = True
185         for icell in xrange(partFixm.getNumberOfCells()):
186             ids5 = iPart.findIdsEqual(icell)
187             areaOfCells = areaPartFixMob[ids5]
188             areaOfCells.abs()
189             if abs(areaOfCells.accumulate()[0] - areaPartFixm[icell]) > 1e-9:
190                 isCheck2OK = False
191                 pass
192             pass
193         print "Check #2? %s" % (str(isCheck2OK))
194
195 Utiliser les informations de l'intersection pour en faire des champs
196 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
197
198 OK pour ``partFixMob``. Nous souhaitons maintenant créer un champ représentant une fonction indicatrice de la zone 
199
200 Maintenant créer un champ aux cellules sur ``partFixMob`` en mettant 0 sur la partie 
201 exclusive ``partFixm`` et 1 sur la partie couverte. Nous créons donc un champ représentant une fonction indicatrice. 
202 Le visualiser en utilisant un fichier VTK (ne pas oublier l'option *Triangulate* de ParaView). ::
203
204         f = ml.MEDCouplingFieldDouble(ml.ON_CELLS,ml.ONE_TIME)
205         m = partFixMob.deepCopy()
206         m.tessellate2D(0.1)
207         f.setMesh(m)
208         arr = ml.DataArrayDouble(partFixMob.getNumberOfCells(),1)
209         arr[iMob.findIdsEqual(-1)] = 0.
210         arr[iMob.findIdsNotEqual(-1)] = 1.
211         f.setArray(arr)
212         f.checkConsistencyLight()
213         f.setName("Zone")
214         ml.MEDCouplingFieldDouble.WriteVTK("Zone.vtu",[f])
215
216 .. image:: images/LocationEx2.jpg
217         :scale: 100
218
219 Plus généralement prendre les zones 0, 1 et 5. Faire un champ aux cellules qui vaut 0 dans la zone exclusivement de ``fixm``,
220 1 dans zone #0, 2 dans la zone #1 et finalement 3 dans la zone #5. ::
221
222         zonesMobm = ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords([mobm[zonesInMobm[0]], mobm[zonesInMobm[1]], mobm[zonesInMobm[5]]])
223         zonesMobm.zipCoords()
224         partFixMob2,iPart2,iMob2 = ml.MEDCouplingUMesh.Intersect2DMeshes(partFixm,zonesMobm,1e-10)
225         partFixMob2.mergeNodes(1e-10)
226         f2 = ml.MEDCouplingFieldDouble(ml.ON_CELLS, ml.ONE_TIME)
227         m2 = partFixMob2.deepCopy()
228         m2.tessellate2D(0.1)
229         f2.setMesh(m2)
230         arr = ml.DataArrayDouble(partFixMob2.getNumberOfCells(),1)
231         arr[iMob2.findIdsEqual(-1)]=0.
232         st = 0
233         end = st + len(zonesInMobm[0])
234         arr[iMob2.findIdsInRange(st,end)] = 1.
235         st += len(zonesInMobm[0]) ; 
236         end = st + len(zonesInMobm[1])
237         arr[iMob2.findIdsInRange(st,end)] = 2.
238         st += len(zonesInMobm[1])
239         end = st + len(zonesInMobm[2])
240         arr[iMob2.findIdsInRange(st,end)] = 3.
241         f2.setArray(arr)
242         f2.checkConsistencyLight()
243         f2.setName("Zone2")
244         ml.MEDCouplingFieldDouble.WriteVTK("Zone2.vtu",[f2])
245
246 Ne pas oublier l'option *Triangulate* de ParaView dans le panneau Display pour bien voir les champs:
247
248 .. image:: images/zonesMobm.jpg
249
250 Solution
251 ~~~~~~~~
252
253 :ref:`python_testmedcouplingloaderex2_solution`