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