doc/tutorial/medcoupling_umesh1_fr.rst

   1
   2 Manipuler les maillages non structurés
   3 --------------------------------------
   4
   5 Les meshes non-structurées sont le type de maillage le plus utilisé. ``MEDCouplingUMesh`` est le nom de la classe en charge
   6 de représenter ces maillages dans MEDCoupling. ``MEDCouplingUMesh`` hérite de la classe ``MEDCouplingPointSet``.
   7 ``MEDCouplingPointSet`` gère toutes les méthodes relatives au coordonnées. ``MEDCouplingUMesh`` a deux attributs en plus de
   8 ceux de ``MEDCouplingPointSet`` permettant de décrire la liste des noeuds contribuants à une cellule (i.e. la *connectivité*).
   9
  10 Objectifs
  11 ~~~~~~~~~
  12
  13 Le but ici est de manipuler des maillages non structurés (en extraire une partie, etc...).
  14 Plusieurs points seront traités dans cet exercice :
  15
  16 * modification des coordonnées d'un maillage
  17 * extraction d'une coupe d'un maillage
  18 * extraire une partie de maillage à partir d'identifiants de cellules
  19 * manipuler les indices, etc ...
  20 * manipulation de la connectivité descendante
  21
  22 .. image:: images/UMesh1.png
  23         :scale: 80
  24
  25 Début de l'implémentation
  26 ~~~~~~~~~~~~~~~~~~~~~~~~~
  27
  28 Importer le module Python ``MEDCoupling``. ::
  29
  30         import MEDCoupling as mc
  31
  32 Construire un maillage. Ce maillage ``mesh3D`` contient artificiellement 2 types de cellules (``mc.NORM_HEXA8`` et ``mc.NORM_POLYHED``)
  33 pour appréhender le mélange de types geometriques.
  34 ``mesh3D`` est un *maillage extrudé* contenant 18 cellules composées de 3 niveaux selon Z, chaque niveau ayant 6 cellules.
  35 Faire un bon gros copier-coller des lignes suivantes pour construire la mesh (l'intérêt de l'exercise vient après) : ::
  36
  37         coords=[0.,0.,0., 1.,1.,0., 1.,1.25,0., 1.,0.,0., 1.,1.5,0., 2.,0.,0., 2.,1.,0., 1.,2.,0., 0.,2.,0., 3.,1.,0.,
  38                 3.,2.,0., 0.,1.,0., 1.,3.,0., 2.,2.,0., 2.,3.,0.,
  39                 0.,0.,1., 1.,1.,1., 1.,1.25,1., 1.,0.,1., 1.,1.5,1., 2.,0.,1., 2.,1.,1., 1.,2.,1., 0.,2.,1., 3.,1.,1.,
  40                 3.,2.,1., 0.,1.,1., 1.,3.,1., 2.,2.,1., 2.,3.,1.,
  41                 0.,0.,2., 1.,1.,2., 1.,1.25,2., 1.,0.,2., 1.,1.5,2., 2.,0.,2., 2.,1.,2., 1.,2.,2., 0.,2.,2., 3.,1.,2.,
  42                 3.,2.,2., 0.,1.,2., 1.,3.,2., 2.,2.,2., 2.,3.,2.,
  43                 0.,0.,3., 1.,1.,3., 1.,1.25,3., 1.,0.,3., 1.,1.5,3., 2.,0.,3., 2.,1.,3., 1.,2.,3., 0.,2.,3., 3.,1.,3.,
  44                 3.,2.,3., 0.,1.,3., 1.,3.,3., 2.,2.,3., 2.,3.,3.]
  45         conn=[0,11,1,3,15,26,16,18,   1,2,4,7,13,6,-1,1,16,21,6,-1,6,21,28,13,-1,13,7,22,28,-1,7,4,19,22,-1,4,2,17,19,-1,2,1,16,17,-1,16,21,28,22,19,17,
  46               1,6,5,3,16,21,20,18,   13,10,9,6,28,25,24,21, 11,8,7,4,2,1,-1,11,26,16,1,-1,1,16,17,2,-1,2,17,19,4,-1,4,19,22,7,-1,7,8,23,22,-1,8,11,26,23,-1,26,16,17,19,22,23,
  47               7,12,14,13,22,27,29,28,  15,26,16,18,30,41,31,33, 16,17,19,22,28,21,-1,16,31,36,21,-1,21,36,43,28,-1,28,22,37,43,-1,22,19,34,37,-1,19,17,32,34,-1,17,16,31,32,-1,31,36,43,37,34,32,
  48               16,21,20,18,31,36,35,33,   28,25,24,21,43,40,39,36, 26,23,22,19,17,16,-1,26,41,31,16,-1,16,31,32,17,-1,17,32,34,19,-1,19,34,37,22,-1,22,23,38,37,-1,23,26,41,38,-1,41,31,32,34,37,38,
  49               22,27,29,28,37,42,44,43, 30,41,31,33,45,56,46,48,  31,32,34,37,43,36,-1,31,46,51,36,-1,36,51,58,43,-1,43,37,52,58,-1,37,34,49,52,-1,34,32,47,49,-1,32,31,46,47,-1,46,51,58,52,49,47,
  50               31,36,35,33,46,51,50,48,  43,40,39,36,58,55,54,51, 41,38,37,34,32,31,-1,41,56,46,31,-1,31,46,47,32,-1,32,47,49,34,-1,34,49,52,37,-1,37,38,53,52,-1,38,41,56,53,-1,56,46,47,49,52,53,
  51               37,42,44,43,52,57,59,58]
  52         mesh3D = mc.MEDCouplingUMesh("mesh3D",3)
  53         mesh3D.allocateCells(18)
  54         mesh3D.insertNextCell(mc.NORM_HEXA8,conn[0:8]); mesh3D.insertNextCell(mc.NORM_POLYHED,conn[8:51]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[51:59]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[59:67]); mesh3D.insertNextCell(mc.NORM_POLYHED,conn[67:110]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[110:118]);
  55         mesh3D.insertNextCell(mc.NORM_HEXA8,conn[118:126]); mesh3D.insertNextCell(mc.NORM_POLYHED,conn[126:169]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[169:177]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[177:185]); mesh3D.insertNextCell(mc.NORM_POLYHED,conn[185:228]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[228:236]);
  56         mesh3D.insertNextCell(mc.NORM_HEXA8,conn[236:244]); mesh3D.insertNextCell(mc.NORM_POLYHED,conn[244:287]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[287:295]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[295:303]); mesh3D.insertNextCell(mc.NORM_POLYHED,conn[303:346]); mesh3D.insertNextCell(mc.NORM_HEXA8,conn[346:354]);
  57         myCoords = mc.DataArrayDouble(coords,60,3)
  58         myCoords.setInfoOnComponents(["X [m]","Y [m]","Z [m]"])
  59         mesh3D.setCoords(myCoords)
  60         mesh3D.orientCorrectlyPolyhedrons()
  61         mesh3D.sortCellsInMEDFileFrmt()
  62         mesh3D.checkConsistencyLight()
  63         renum = mc.DataArrayInt(60); renum[:15]=range(15,30) ; renum[15:30]=range(15) ; renum[30:45]=range(45,60) ; renum[45:]=range(30,45)
  64         mesh3D.renumberNodes(renum,60)
  65
  66 Convertir les unités
  67 ~~~~~~~~~~~~~~~~~~~~
  68
  69 On convertit ici les coordonnées de mètres en centimètres.
  70 Cela paraît idiot mais c'est un très grand classique du couplage ... ::
  71
  72         mesh3D.getCoords()[:] *= 100.
  73         mesh3D.getCoords().setInfoOnComponents(["X [cm]","Y [cm]","Z [cm]"])
  74
  75 .. note:: Il est important de mettre à jour les informations sur les composantes des coordonnées (les unités) pour éviter toute ambiguïté.
  76         INTERP_KERNEL library inclut un évaluateur d'unité.
  77
  78 .. note:: Noter l'astuce sur la première ligne ``[:]`` afin de récupérer la version inscriptible des coordonnées
  79         (et non une copie temporaire)
  80
  81 Trouver les différents niveaux
  82 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  83
  84 Le maillage est extrudé, il est donc très régulier, et aligné sur les axes Ox, Oy et Oz (cf figure).
  85 On veut connaître quelles
  86 sont les côtes Z des différentes couches de cubes.
  87 Extraire les différents niveaux en Z dans ``mesh3D``, rangés de manière croissante.
  88 Utiliser la méthode ``DataArrayDouble.getDifferentValues()`` and ``DataArrayDouble.sort()``. ::
  89
  90         zLev = mesh3D.getCoords()[:,2]
  91         zLev = zLev.getDifferentValues(1e-12)
  92         zLev.sort()     # In-place sort
  93
  94 Extraire des identifiants de cellules
  95 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  96
  97 Extraire les 6 identifiants des cellules de la seconde rangée suivant Oz.
  98 Il y a 3 possibilités pour faire cela. Nous allons les voir du plus simple au plus complexe.
  99
 100 * En utilisant ``buildSlice3D()`` :
 101         Méthode très simple mais gourmande en CPU. Pour trouver la solution il suffit de définir un plan dont le vecteur normal est ``[0.,0.,1.]``
 102         et passant par le point ``[0., 0., (zLev[1]+zLev[2])/2]``.
 103         La méthode retourne deux choses : le maillage de coupe ``tmp`` (un maillage de mesh-dimension 2, mais de dimension spatiale
 104         3) et pour chaque cellule 3D surfacique de ``tmp``, l'identifiant de la cellule 3D (=un volume) coupée dans le
 105         maillage de départ  ::
 106
 107                 tmp, cellIdsSol1 = mesh3D.buildSlice3D([0.,0.,(zLev[1]+zLev[2])/2], [0.,0.,1.], 1e-12)
 108
 109 * En utilisant les barycentres des cellules de ``mesh3D`` :
 110         L'utilisation des barycentres est une technique classique pour identifier un ensemble de cellules répondant à certains
 111         critères géométriques.
 112         Il s'agit d'abord de calculer les barycentres des cellules 3D de ``mesh3D`` (méthode
 113         ``MEDCouplingUMesh.computeCellCenterOfMass()``).
 114         (*Note*: le nom -- un peu trop long -- de cette méthode hérite du passé. Le "AndOwner" indique le fait qu'en C++
 115         l'appelant est responsable de la désallocation de l'objet retourné : il prend l'*ownership* du résultat).
 116
 117         Ensuite sélectionner la composante #2 des barycentres des cellules et mettre le résultat dans ``baryZ``.
 118         Ensuite il suffit de selectionner dans ``baryZ`` les tuples qui sont dans l'intervalle ``[zLev[1], zLev[2]]``.
 119         Les identifiants de ces tuples (i.e. leur index dans ``baryZ``) est directement un identifiant de cellule
 120         car ``computeCellCenterOfMass()`` retourne un tableau indéxé par les numéros de cellule.::
 121
 122                 bary = mesh3D.computeCellCenterOfMass()
 123                 baryZ = bary[:,2]
 124                 cellIdsSol2 = baryZ.findIdsInRange(zLev[1], zLev[2])
 125
 126 * En utilisant ``MEDCouplingMappedExtrudedMesh`` :
 127         C'est la méthode exclusivement basée sur la connectivité nodale pour déduire l'extrusion. Les coordonnées sont ici ignorées.
 128         Pour construire un ``MEDCouplingMappedExtrudedMesh`` deux objets sont requis. Le maillage non-structuré 3D
 129         représentant en fait un maillage *extrudé*, et un maillage non structuré 3D surfacique (mesh-dim 2)
 130         reposant sur les mêmes coordonnéees, à partir duquel l'extrusion sera calculée.
 131         Commencer par construire le maillage 3D surfacique. Pour ce faire il suffit de repérer les noeuds appartenant
 132         à 1e-10 près de plan de vecteur normal ``[0.,0.,1.]`` et passant
 133         par ``[0.,0.,zLev[0]]`` (``MEDCouplingUMesh.findNodesOnPlane()``). Ensuite appeler ``MEDCouplingUMesh.buildFacePartOfMySelfNode()``
 134         pour construire ``mesh2D`` (lire la doc de la fonction). ::
 135
 136                 nodeIds = mesh3D.findNodesOnPlane([0., 0., zLev[0]], [0.,0.,1.], 1e-10)
 137                 mesh2D = mesh3D.buildFacePartOfMySelfNode(nodeIds, True)
 138
 139
 140         Il est alors possible de construire un maillage extrudé ``extMesh`` à partir de ``mesh3D`` et de ``mesh2D``.
 141         Un maillage extrudé se construit en *reconnaissant* un maillage non structuré comme étant l'extrusion d'un maillage
 142         de dimension ``n-1`` (avec ``n`` la dimension initiale de ``mesh3D``, ici 3). Si cela n'est pas le cas, la construction
 143         plante. Le maillage 2D est forcément en haut ou en bas du 3D volumique, et le dernier entier spécifie la cellule à partir
 144         de laquelle le fil de fer 1D guidant l'extrusion sera construit : ::
 145
 146                 extMesh = mc.MEDCouplingMappedExtrudedMesh(mesh3D, mesh2D, 0)
 147
 148         On a alors la garantie que, dans ``extMesh``,  les cellules sont ordonnées par niveau Z croissant.
 149         Il suffit de récupérer le 2ème niveau (``MEDCouplingMappedExtrudedMesh.getMesh3DIds()``). ::
 150
 151                 n_cells = mesh2D.getNumberOfCells()
 152                 cellIdsSol3 = extMesh.getMesh3DIds()[n_cells:2*n_cells]
 153
 154 On vérifie alors que les 3 solutions sont les mêmes : ::
 155
 156         print cellIdsSol1.getValues()
 157         print cellIdsSol2.getValues()
 158         print cellIdsSol3.getValues()
 159
 160
 161 Extraire une sous partie d'un maillage 3D
 162 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 163
 164 Utiliser les identifiants de cellules ``cellIdsSol2`` obtenus précédemment pour extraire une sous-partie de ``mesh3D``,
 165 c'est-à-dire un maillage avec un sous-ensemble des cellules de ``mesh3D``. ::
 166
 167         mesh3DPart = mesh3D[cellIdsSol2]
 168
 169 .. note:: En C++ la méthode sous-jacente invoquée (et par ailleurs aussi disponible en Python) s'appelle
 170         ``mesh3DPart = mesh3D.buildPartOfMySelf(cellIdsSol2,True)``
 171
 172 .. note:: Le type géométrique ne rentre pas du tout en compte ici. L'instruction précédente prend les cellules
 173         dans l'ordre où elles sont disponibles dans le maillage initial.
 174
 175 L'objet ``mesh3DPart`` contient ``len(cellIdsSol2)`` cellules désormais. La cellule #0 de ``mesh3DPart`` correspond à la cellule avec l'identifiant ``cellIdsSol2[0]`` de ``mesh3D``, et ainsi de suite. Ainsi ``cellIdsSol2`` peut être vu comme un
 176 tableau new-2-old.
 177
 178 A ce point, ``mesh3DPart`` repose sur une copie du tableau de coordonnées de ``mesh3D``, c'est-à-dire  60 nodes.
 179 Seuls 30 sont effectivement utilisés.
 180 Pour retirer les noeuds orphelins de ``mesh3DPart`` invoquer simplement ``MEDCouplingUMesh.zipCoords()``. ::
 181
 182         mesh3DPart.zipCoords()
 183
 184 Maintenant, ``mesh3DPart`` repose sur 30 nodes et possède 6 cellules. Pour être prêt aux I/O MED-fichier, il est
 185 alors important de voir si ``mesh3DPart`` est bien ordonné, c'est-à-dire si ses cellules sont bien rangées par type géométrique.
 186 On commence par inspecter l'état actuel : ::
 187
 188         print mesh3DPart.advancedRepr()
 189
 190 La fonction suivante fait le même travail : ::
 191
 192         print mesh3DPart.checkConsecutiveCellTypesAndOrder([mc.NORM_HEXA8, mc.NORM_POLYHED])
 193
 194 Ou bien : ::
 195
 196         print mesh3DPart.checkConsecutiveCellTypes()
 197
 198 On voit que ``mesh3DPart`` contient 6 cellules, quatre HEXA8 puis deux POLYHED. Les cellules sont bien
 199 groupées par type géométrique. Si ce n'était pas le cas, on aurait pu invoquer ``MEDCouplingUMesh.sortCellsInMEDFileFrmt()``.
 200
 201
 202 Extraire des cellules alignées sur une ligne 3D
 203 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 204
 205 On souhaite extraire de ``mesh3D`` les 3 cellules dont les barycentres sont le long de la ligne portée par
 206 ``v = [0.,0.,1.]`` et passant par ``pt = [250.,150.,0.]``.
 207 Il y a deux solutions.
 208
 209 * les barycentres de ``mesh3D``  : même principe qu'au-dessus. ::
 210
 211         baryXY = bary[:,[0,1]]
 212         baryXY -= [250.,150.]
 213         magn = baryXY.magnitude()
 214         cellIds2Sol1 = magn.findIdsInRange(0.,1e-12)
 215
 216 * utiliser le maillage extrudé ``extMesh`` : partant de l'unique cellule dans ``mesh2D`` dont le centre est
 217   en ``[250.,150.,0.]``, la méthdode ``MEDCouplingMappedExtrudedMesh.getMesh3DIds()`` retourne les identifiants de
 218   cellules rangée par rangée. ::
 219
 220         bary2 = mesh2D.computeCellCenterOfMass()[:,[0,1]]
 221         bary2 -= [250.,150.]
 222         magn = bary2.magnitude()
 223         ids = magn.findIdsInRange(0.,1e-12)
 224         idStart = int(ids) # ids is assumed to contain only one value, if not an exception is thrown
 225         ze_range = range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells())
 226         cellIds2Sol2 = extMesh.getMesh3DIds()[ze_range]
 227
 228 Maintenant on construit cette sous partie de ``mesh3D`` en utilisant ``cellIds2Sol1`` ou ``cellIds2Sol2``: ::
 229
 230         mesh3DSlice2 = mesh3D[cellIds2Sol1]
 231         mesh3DSlice2.zipCoords()
 232
 233 Duplication, translation et aggrégation de maillages
 234 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 235
 236 Cette partie de l'exercice est intéressante pour construire des maillages complexes, ou pour aggréger des parties
 237 de maillages venant de différents processeurs.
 238
 239 On cherche ici à dupliquer ``mesh3DSlice2``, le translater et l'aggréger avec l'original.
 240
 241 Effectuer une copie complète de ``mesh3DSlice2`` (aussi appelée *deep copy*) sous le nom ``mesh3DSlice2bis``.
 242 Sur cette copie effectuer une translation de ``v=[0.,1000.,0.]``.
 243 Puis aggréger ``mesh3DSlice2`` avec sa copie translatée ``mesh3DSlice2bis``, en utilisant ``MEDCouplingUMesh.MergeUMeshes()``. ::
 244
 245         mesh3DSlice2bis = mesh3DSlice2.deepCopy()
 246         mesh3DSlice2bis.translate([0.,1000.,0.])
 247         mesh3DSlice2All = mc.MEDCouplingUMesh.MergeUMeshes([mesh3DSlice2,mesh3DSlice2bis])
 248         mesh3DSlice2All.writeVTK("mesh3DSlice2All.vtu")
 249
 250 .. note:: Pour information pour merger deux (ou plus) maillages non structurés, il faut invoquer ``MEDCouplingUMesh.MergeUMeshes()``
 251         puis ``MEDCouplingUMesh.mergeNodes()`` sur le résultat, et enfin ``MEDCouplingUMesh.zipConnectivity()``.
 252
 253 .. _exo-umesh-desc-connec:
 254
 255 Connectivité descendante
 256 ~~~~~~~~~~~~~~~~~~~~~~~~
 257
 258 Le but ici est de présenter la notion de *connectivité descendante* (*descending connectivity*).
 259
 260 La connectivité descendante représente les éléments de dimension ``n-1``
 261 constituant chacune des cellules de dimension ``n`` (avec donc ``n`` la dimension du maillage, *mesh-dim*). Par exemple, pour un
 262 maillage de dimension 3 (les cellules sont des *volumes* 3D), cela donne l'ensemble des faces (des *surfaces* 2D) bordant
 263 ces volumes.
 264
 265 A titre d'exemple, on se propose dans notre cas de récupérer les faces *internes* du maillage ``mesh3D``.
 266 Pour cela il est nécessaire de construire le maillage
 267 descendant de ``mesh3D`` (stocké dans ``mesh3DSurf``) c'est-à-dire
 268 le maillage de mesh-dimension 2 (soit ``mesh3D.getMeshDimension()-1``) constitué
 269 des *faces* bordant chacune des cellules (ici des *volumes* 3D) de ``mesh3D``.
 270 La méthode ``MEDCoupling.buildDescendingConnectivity()`` calcule ce maillage, et retourne en même temps des tableaux
 271 de correspondance. Ces tableaux font le lien entre les identifiants des cellules de ``mesh3D``
 272 vers les identifiants de cellules de ``mesh3DSurf``, et vice-et-versa.
 273
 274 Une face de ``mesh3DSurf`` est dite interne, si et seulement si, elle est partagée par plus d'une cellule 3D de ``mesh3D``.
 275 Les 3ème et 4ème paramètres de sortie de la fonction donnent le lien
 276 entre une face et ses cellules *parentes* (i.e. le ou les volumes qu'elle délimite).
 277 Ce lien est exprimé au format *indirect index* vu dans le premier exercice :ref:`indirect-index-exo`. ::
 278
 279         mesh3DSurf, desc, descIndx, revDesc, revDescIndx = mesh3D.buildDescendingConnectivity()
 280         numberOf3DCellSharing = revDescIndx.deltaShiftIndex()
 281         cellIds = numberOf3DCellSharing.findIdsNotEqual(1)
 282         mesh3DSurfInside = mesh3DSurf[cellIds]
 283         mesh3DSurfInside.writeVTK("mesh3DSurfInside.vtu")
 284
 285 Ce genre de manipulation est très utile pour accéder au voisinage d'une ou plusieurs cellules d'un maillage non-structuré.
 286
 287 .. image:: images/mesh3DSurfInside.jpg
 288
 289 Solution
 290 ~~~~~~~~
 291
 292 :ref:`python_testMEDCouplingumesh1_solution`