2 MEDCoupling / MEDLoader - Exemple complet 2 - RJH
3 -------------------------------------------------
5 Ici nous partons de deux fichiers MED très librement inspirés du réacteur expérimental
6 RJH (Réacteur Jules Horowitz).
8 Le premier, :download:`Fixe.med <data/Fixe.med>`, représente la partie statique du réacteur
9 sans dispositif expérimental.
11 .. image:: images/fixm.jpg
14 Le deuxième, :download:`Mobile.med <data/Mobile.med>`, représente la partie mobile.
16 .. image:: images/mobm.jpg
23 Le but ici est d'utiliser MEDCoupling pour:
25 * intersecter ces deux maillages,
27 * et ainsi localiser les zones de recouvrement
30 Début de l'implémentation
31 ~~~~~~~~~~~~~~~~~~~~~~~~~
33 Cet exercice repose comme tous les autres sur le language de script Python. On charge
34 le module Python ``medcoupling``.::
36 import medcoupling as mc
38 Lire et réparer le maillage statique "Fixe.med"
39 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
41 Avec l'API avancée lire tout le fichier "Fixe.med" et appeler ``fixm``
42 l'objet de type ``MEDCouplingUMesh`` représentant le maillage statique. ::
44 fixe = mc.MEDFileMesh.New("Fixe.med")
45 fixm = fixe.getMeshAtLevel(0)
47 Pour ce qui suit il faut absolument que deux cellules se touchant partagent les mêmes edges. Pour ce faire, comme on est
48 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.
49 Fusionner le noeuds distants de moins de 1e-10 et regarder l'impact sur le nombre de noeuds de ``fixm``. ::
51 print("Nb of nodes in the file : %i " % (fixm.getNumberOfNodes()))
52 fixm.mergeNodes(1e-10)
53 print("Nb of non duplicated nodes : %i" % (fixm.getNumberOfNodes()))
55 Même traitement pour ``Mobile.med``, le lire avec l'API avancée de MEDLoader (appeler ``mobm`` l'instance du maillage)
56 et le réparer en supprimant les noeuds dupliqués. ::
58 mobile = mc.MEDFileMesh.New("Mobile.med")
59 mobm = mobile.getMeshAtLevel(0)
60 mobm.mergeNodes(1e-10)
62 Le maillage du RJH étant plus général que des ``TRI6`` et des ``QUAD8``, on a besoin
63 de stocker ces cellules avec un type géométrique à ``QPOLYG`` (Quadratic Polygon) qui représente un polygone *quadratique*
64 (le terme n'est pas très heureux, encore des raisons historiques, ...), c'est-à-dire un polygone avec un nombre arbitraire
65 de côtés, et potentiellement des côtés en forme d'arcs de cercles plutôt que de segments de droites.
66 Ce type géométrique ``NORM_QPOLYG`` est dans MEDCoupling/MEDLoader et aussi dans MED fichier.
68 Nous voudrions visualiser ces deux maillages dans PARAVIS/ParaView, mais nous rencontrons ici deux soucis:
70 * les polygones non-convexes sont, par défaut, mal représentés par VTK en mode *Surface*.
71 Il faut sélectionner l'option avancée "Triangulate" dans le panneau Display de PARAVIS/ParaView pour avoir un rendu correct.
72 * les arcs de cercles ne sont pas correctement supportés par ParaView. Il faut les *tesséliser*, c'est-à-dire les transformer
73 en plusieurs petits segments de droite. La méthode ``MEDCouplingUMesh.tessellate2D()`` fait ce travail, mais modifie
74 le maillage. Nous faisons donc une copie préalable. Nous passons en paramètre la finesse de découpage (0.1 est suffisant
75 -- angle en radian). Attention donc à ne pas modifer ni ``fixm`` ni ``mobm`` ! ::
77 fixm2 = fixm.deepCopy() # tessellate2D() modifies the current mesh
78 fixm2.tessellate2D(0.1)
79 fixm2.writeVTK("fixm2.vtu")
80 mobm2 = mobm.deepCopy()
81 mobm2.tessellate2D(0.1)
82 mobm2.writeVTK("mobm2.vtu")
84 Faire une petite méthode ``displayVTK()``, faisant le travail qui nous servira souvent après. ::
86 def displayVTK(m,fname):
92 Faire des réductions et des repérages de zones
93 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
95 Le maillage ``mobm`` est en 6 parties distinctes (voir l'image au dessus). On ne veut récupérer que la première partie.
96 Utiliser ``MEDCouplingUMesh.partitionBySpreadZone()`` pour partitionner en zones ``mobm`` et ne prendre que la première zone.
97 Appeler cette nouvelle instance ``zone1Mobm`` et lui retirer tous les noeuds orphelins (``MEDCouplingUMesh.zipCoords()``)
100 zonesInMobm = mobm.partitionBySpreadZone()
101 print("Nb of zones in mobm : %i" % (len(zonesInMobm)))
102 zone1Mobm = mobm[zonesInMobm[0]]
103 zone1Mobm.zipCoords()
104 displayVTK(zone1Mobm, "zone1Mobm.vtu")
106 .. image:: images/zone1Mobm.jpg
109 Nous allons désormais travailler autour de ``zone1Mobm``. Nous allons réduire la zone de travail de ``fixm`` autour de ``zone1Mobm``.
110 Pour ce faire, réduire ``fixm`` en ne prenant que les cellules dans la boîte englobante
111 de ``zone1Mobm`` (``MEDCouplingUMesh.getBoundingBox()`` et ``MEDCouplingUMesh.getCellsInBoundingBox()``).
112 Appeler ce nouvel objet ``partFixm``, lui retirer ses noeuds orphelins et l'afficher. ::
114 ids2 = fixm.getCellsInBoundingBox(zone1Mobm.getBoundingBox(),1e-10)
115 partFixm = fixm[ids2]
117 displayVTK(partFixm,"partFixm.vtu")
119 .. image:: images/partFixmAndzone1Mobm.jpg
121 Intersecter géométriquement deux maillages
122 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
124 C'est le coeur de l'exercice. Nous allons intersecter géométriquement ``partFixm`` et ``zone1Mobm``. Cela revient à
125 partitionner à minima ``partFixm`` en cellules appartenant
126 soit complètement à ``partFixm`` soit à ``partFixm`` et ``zone1Mobm``. Invoquer la méthode statique
127 ``MEDCouplingUMesh.Intersect2DMeshes()``, avec ``partFixm`` et ``zone1Mobm`` et mettre une précision
128 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
129 ici ``partFixMob``, ``iPart`` et ``iMob`` dans cet ordre.
131 Sur ``partFixMob`` merger les noeuds à 1e-10 près. ::
133 partFixMob, iPart, iMob = mc.MEDCouplingUMesh.Intersect2DMeshes(partFixm,zone1Mobm,1e-10)
134 partFixMob.mergeNodes(1e-10)
136 Récupérer et afficher la partie de ``partFixm`` qui n'est pas dans ``zone1Mobm``. Appeler ce maillage ``partFixmWithoutZone1Mobm``. ::
138 ids3 = iMob.findIdsEqual(-1)
139 partFixmWithoutZone1Mobm = partFixMob[ids3]
140 displayVTK(partFixmWithoutZone1Mobm,"partFixmWithoutZone1Mobm.vtu")
142 .. image:: images/partFixmWithoutZone1Mobm.jpg
145 Maintenant, on va vérifier la qualité du résultat retourné par ``MEDCouplingUMesh.Intersect2DMeshes()``.
146 Pour ce faire on va passer 3 tests:
148 * **Check #0** la somme des aires des cellules de ``partFixm`` et égale à celle de ``partFixMob``
149 * **Check #1** la somme des aires des cellules de ``zone1Mobm`` et égale à la somme des cells de ``partFixMob``
150 dont l'id dans ``iMob`` est different de -1
151 * **Check #2** pour chaque cellule de ``partFixm``, son aire est égale à la somme des aires des cellules de ``partFixMob``
153 L'aire est une valeur algébrique. Donc attention cette verification ne peut se faire que si les cellules
154 sont toutes bien orientées ou à minima toutes orientées de la même manière.
155 Pour ce faire, regardons les aires des 38 cellules de ``partFixm`` (nom de variable : ``areaPartFixm``). ::
157 areaPartFixm = partFixm.getMeasureField(True).getArray()
158 print(areaPartFixm.getValues())
160 On voit que toutes les valeurs sont négatives. *Bilan*: ce fichier MED ne respecte pas la convention MED fichier !
161 ``partFixm`` étant mal orienté, et ``MEDCouplingUMesh.Intersect2DMeshes()`` conservant l'orientation,
162 ``partFixMob`` est lui aussi mal orienté.
163 Bref, on va faire les comparaisons sur des tableaux de valeurs absolues. Vérifier alors **Check #0**. ::
165 areaPartFixm = partFixm.getMeasureField(isAbs=False).getArray() # prise en compte de l'orientation des mailles
167 areaPartFixMob = partFixMob.getMeasureField(isAbs=False).getArray()
169 val1=areaPartFixm.accumulate()[0]
170 val2=areaPartFixMob.accumulate()[0]
171 print("Check #0 %lf == %lf with precision 1e-8? %s" % (val1,val2,str(abs(val1-val2)<1e-8)))
173 On peut passer au **Check #1**. L'esprit est le même que le **Check #0**. ::
175 areaZone1Mobm = zone1Mobm.getMeasureField(isAbs=False).getArray()
177 val3 = areaZone1Mobm.accumulate()[0]
178 ids4 = iMob.findIdsNotEqual(-1)
179 areaPartFixMob2 = areaPartFixMob[ids4]
180 val4 = areaPartFixMob2.accumulate()[0]
181 print("Check #1 %lf == %lf with precision 1e-8 ? %s" % (val3,val4,str(abs(val3-val4)<1e-8)))
183 Puis le **Check #2**. ::
186 for icell in list(range(partFixm.getNumberOfCells())):
187 ids5 = iPart.findIdsEqual(icell)
188 areaOfCells = areaPartFixMob[ids5]
190 if abs(areaOfCells.accumulate()[0] - areaPartFixm[icell]) > 1e-9:
194 print("Check #2? %s" % (str(isCheck2OK)))
196 Utiliser les informations de l'intersection pour en faire des champs
197 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
199 OK pour ``partFixMob``. Nous souhaitons maintenant créer un champ représentant une fonction indicatrice de la zone
201 Maintenant créer un champ aux cellules sur ``partFixMob`` en mettant 0 sur la partie
202 exclusive ``partFixm`` et 1 sur la partie couverte. Nous créons donc un champ représentant une fonction indicatrice.
203 Le visualiser en utilisant un fichier VTK (ne pas oublier l'option *Triangulate* de ParaView). ::
205 f = mc.MEDCouplingFieldDouble(mc.ON_CELLS,mc.ONE_TIME)
206 m = partFixMob.deepCopy()
209 arr = mc.DataArrayDouble(partFixMob.getNumberOfCells(),1)
210 arr[iMob.findIdsEqual(-1)] = 0.
211 arr[iMob.findIdsNotEqual(-1)] = 1.
213 f.checkConsistencyLight()
215 mc.MEDCouplingFieldDouble.WriteVTK("Zone.vtu",[f])
217 .. image:: images/LocationEx2.jpg
220 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``,
221 1 dans zone #0, 2 dans la zone #1 et finalement 3 dans la zone #5. ::
223 zonesMobm = mc.MEDCouplingUMesh.MergeUMeshesOnSameCoords([mobm[zonesInMobm[0]], mobm[zonesInMobm[1]], mobm[zonesInMobm[5]]])
224 zonesMobm.zipCoords()
225 partFixMob2,iPart2,iMob2 = mc.MEDCouplingUMesh.Intersect2DMeshes(partFixm,zonesMobm,1e-10)
226 partFixMob2.mergeNodes(1e-10)
227 f2 = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
228 m2 = partFixMob2.deepCopy()
231 arr = mc.DataArrayDouble(partFixMob2.getNumberOfCells(),1)
232 arr[iMob2.findIdsEqual(-1)]=0.
234 end = st + len(zonesInMobm[0])
235 arr[iMob2.findIdsInRange(st,end)] = 1.
236 st += len(zonesInMobm[0]) ;
237 end = st + len(zonesInMobm[1])
238 arr[iMob2.findIdsInRange(st,end)] = 2.
239 st += len(zonesInMobm[1])
240 end = st + len(zonesInMobm[2])
241 arr[iMob2.findIdsInRange(st,end)] = 3.
243 f2.checkConsistencyLight()
245 mc.MEDCouplingFieldDouble.WriteVTK("Zone2.vtu",[f2])
247 Ne pas oublier l'option *Triangulate* de ParaView dans le panneau Display pour bien voir les champs:
249 .. image:: images/zonesMobm.jpg
254 :ref:`python_testmedcouplingloaderex2_solution`