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 Cet exercice repose comme tous les autres sur le language de script Python. On charge
  34 le module Python ``medcoupling``.::
  35
  36     import medcoupling as mc
  37
  38 Lire et réparer le maillage statique "Fixe.med"
  39 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  40
  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. ::
  43
  44         fixe = mc.MEDFileMesh.New("Fixe.med")
  45         fixm = fixe.getMeshAtLevel(0)
  46
  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``. ::
  50
  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()))
  54
  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. ::
  57
  58         mobile = mc.MEDFileMesh.New("Mobile.med")
  59         mobm = mobile.getMeshAtLevel(0)
  60         mobm.mergeNodes(1e-10)
  61
  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.
  67
  68 Nous voudrions visualiser ces deux maillages dans PARAVIS/ParaView, mais nous rencontrons ici deux soucis:
  69
  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`` ! ::
  76
  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")
  83
  84 Faire une petite méthode ``displayVTK()``, faisant le travail qui nous servira souvent après. ::
  85
  86         def displayVTK(m,fname):
  87                 tmp = m.deepCopy()
  88                 tmp.tessellate2D(0.1)
  89                 tmp.writeVTK(fname)
  90                 return
  91
  92 Faire des réductions et des repérages de zones
  93 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  94
  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()``)
  98 Enfin l'afficher : ::
  99
 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")
 105
 106 .. image:: images/zone1Mobm.jpg
 107         :scale: 70
 108
 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. ::
 113
 114         ids2 = fixm.getCellsInBoundingBox(zone1Mobm.getBoundingBox(),1e-10)
 115         partFixm = fixm[ids2]
 116         partFixm.zipCoords()
 117         displayVTK(partFixm,"partFixm.vtu")
 118
 119 .. image:: images/partFixmAndzone1Mobm.jpg
 120
 121 Intersecter géométriquement deux maillages
 122 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 123
 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.
 130
 131 Sur ``partFixMob`` merger les noeuds à 1e-10 près. ::
 132
 133         partFixMob, iPart, iMob = mc.MEDCouplingUMesh.Intersect2DMeshes(partFixm,zone1Mobm,1e-10)
 134         partFixMob.mergeNodes(1e-10)
 135
 136 Récupérer et afficher la partie de ``partFixm`` qui n'est pas dans ``zone1Mobm``. Appeler ce maillage ``partFixmWithoutZone1Mobm``. ::
 137
 138         ids3 = iMob.findIdsEqual(-1)
 139         partFixmWithoutZone1Mobm = partFixMob[ids3]
 140         displayVTK(partFixmWithoutZone1Mobm,"partFixmWithoutZone1Mobm.vtu")
 141
 142 .. image:: images/partFixmWithoutZone1Mobm.jpg
 143         :scale: 70
 144
 145 Maintenant, on va vérifier la qualité du résultat retourné par ``MEDCouplingUMesh.Intersect2DMeshes()``.
 146 Pour ce faire on va passer 3 tests:
 147
 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``
 152
 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``). ::
 156
 157         areaPartFixm = partFixm.getMeasureField(True).getArray()
 158         print(areaPartFixm.getValues())
 159
 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**. ::
 164
 165         areaPartFixm = partFixm.getMeasureField(isAbs=False).getArray() # prise en compte de l'orientation des mailles
 166         areaPartFixm.abs()
 167         areaPartFixMob = partFixMob.getMeasureField(isAbs=False).getArray()
 168         areaPartFixMob.abs()
 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)))
 172
 173 On peut passer au **Check #1**. L'esprit est le même que le **Check #0**. ::
 174
 175         areaZone1Mobm = zone1Mobm.getMeasureField(isAbs=False).getArray()
 176         areaZone1Mobm.abs()
 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)))
 182
 183 Puis le **Check #2**. ::
 184
 185         isCheck2OK = True
 186         for icell in list(range(partFixm.getNumberOfCells())):
 187             ids5 = iPart.findIdsEqual(icell)
 188             areaOfCells = areaPartFixMob[ids5]
 189             areaOfCells.abs()
 190             if abs(areaOfCells.accumulate()[0] - areaPartFixm[icell]) > 1e-9:
 191                 isCheck2OK = False
 192                 pass
 193             pass
 194         print("Check #2? %s" % (str(isCheck2OK)))
 195
 196 Utiliser les informations de l'intersection pour en faire des champs
 197 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 198
 199 OK pour ``partFixMob``. Nous souhaitons maintenant créer un champ représentant une fonction indicatrice de la zone
 200
 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). ::
 204
 205         f = mc.MEDCouplingFieldDouble(mc.ON_CELLS,mc.ONE_TIME)
 206         m = partFixMob.deepCopy()
 207         m.tessellate2D(0.1)
 208         f.setMesh(m)
 209         arr = mc.DataArrayDouble(partFixMob.getNumberOfCells(),1)
 210         arr[iMob.findIdsEqual(-1)] = 0.
 211         arr[iMob.findIdsNotEqual(-1)] = 1.
 212         f.setArray(arr)
 213         f.checkConsistencyLight()
 214         f.setName("Zone")
 215         mc.MEDCouplingFieldDouble.WriteVTK("Zone.vtu",[f])
 216
 217 .. image:: images/LocationEx2.jpg
 218         :scale: 100
 219
 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. ::
 222
 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()
 229         m2.tessellate2D(0.1)
 230         f2.setMesh(m2)
 231         arr = mc.DataArrayDouble(partFixMob2.getNumberOfCells(),1)
 232         arr[iMob2.findIdsEqual(-1)]=0.
 233         st = 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.
 242         f2.setArray(arr)
 243         f2.checkConsistencyLight()
 244         f2.setName("Zone2")
 245         mc.MEDCouplingFieldDouble.WriteVTK("Zone2.vtu",[f2])
 246
 247 Ne pas oublier l'option *Triangulate* de ParaView dans le panneau Display pour bien voir les champs:
 248
 249 .. image:: images/zonesMobm.jpg
 250
 251 Solution
 252 ~~~~~~~~
 253
 254 :ref:`python_testmedcouplingloaderex2_solution`