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.deepCpy()        # tessellate2D() modifies the current mesh
  77         fixm2.tessellate2D(0.1)
  78         fixm2.writeVTK("fixm2.vtu")
  79         mobm2 = mobm.deepCpy()
  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.deepCpy()
  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.getIdsEqual(-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.getIdsNotEqual(-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.getIdsEqual(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.deepCpy()
 206         m.tessellate2D(0.1)
 207         f.setMesh(m)
 208         arr = ml.DataArrayDouble(partFixMob.getNumberOfCells(),1)
 209         arr[iMob.getIdsEqual(-1)] = 0.
 210         arr[iMob.getIdsNotEqual(-1)] = 1.
 211         f.setArray(arr)
 212         f.checkCoherency()
 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.deepCpy()
 228         m2.tessellate2D(0.1)
 229         f2.setMesh(m2)
 230         arr = ml.DataArrayDouble(partFixMob2.getNumberOfCells(),1)
 231         arr[iMob2.getIdsEqual(-1)]=0.
 232         st = 0
 233         end = st + len(zonesInMobm[0])
 234         arr[iMob2.getIdsInRange(st,end)] = 1.
 235         st += len(zonesInMobm[0]) ;
 236         end = st + len(zonesInMobm[1])
 237         arr[iMob2.getIdsInRange(st,end)] = 2.
 238         st += len(zonesInMobm[1])
 239         end = st + len(zonesInMobm[2])
 240         arr[iMob2.getIdsInRange(st,end)] = 3.
 241         f2.setArray(arr)
 242         f2.checkCoherency()
 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`