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