]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/medcoupling_umesh1_fr.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / 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         
115         Ensuite sélectionner la composante #2 des barycentres des cellules et mettre le résultat dans ``baryZ``.
116         Ensuite il suffit de selectionner dans ``baryZ`` les tuples qui sont dans l'intervalle ``[zLev[1], zLev[2]]``. 
117         Les identifiants de ces tuples (i.e. leur index dans ``baryZ``) est directement un identifiant de cellule
118         car ``computeCellCenterOfMass()`` retourne un tableau indéxé par les numéros de cellule.::
119         
120                 bary = mesh3D.computeCellCenterOfMass()
121                 baryZ = bary[:,2]
122                 cellIdsSol2 = baryZ.findIdsInRange(zLev[1], zLev[2])
123
124 * En utilisant ``MEDCouplingMappedExtrudedMesh`` :
125         C'est la méthode exclusivement basée sur la connectivité nodale pour déduire l'extrusion. Les coordonnées sont ici ignorées.
126         Pour construire un ``MEDCouplingMappedExtrudedMesh`` deux objets sont requis. Le maillage non-structuré 3D  
127         représentant en fait un maillage *extrudé*, et un maillage non structuré 3D surfacique (mesh-dim 2) 
128         reposant sur les mêmes coordonnéees, à partir duquel l'extrusion sera calculée.
129         Commencer par construire le maillage 3D surfacique. Pour ce faire il suffit de repérer les noeuds appartenant 
130         à 1e-10 près de plan de vecteur normal ``[0.,0.,1.]`` et passant
131         par ``[0.,0.,zLev[0]]`` (``MEDCouplingUMesh.findNodesOnPlane()``). Ensuite appeler ``MEDCouplingUMesh.buildFacePartOfMySelfNode()`` 
132         pour construire ``mesh2D`` (lire la doc de la fonction). ::
133         
134                 nodeIds = mesh3D.findNodesOnPlane([0., 0., zLev[0]], [0.,0.,1.], 1e-10)
135                 mesh2D = mesh3D.buildFacePartOfMySelfNode(nodeIds, True)
136                 
137
138         Il est alors possible de construire un maillage extrudé ``extMesh`` à partir de ``mesh3D`` et de ``mesh2D``. 
139         Un maillage extrudé se construit en *reconnaissant* un maillage non structuré comme étant l'extrusion d'un maillage
140         de dimension ``n-1`` (avec ``n`` la dimension initiale de ``mesh3D``, ici 3). Si cela n'est pas le cas, la construction
141         plante. Le maillage 2D est forcément en haut ou en bas du 3D volumique, et le dernier entier spécifie la cellule à partir
142         de laquelle le fil de fer 1D guidant l'extrusion sera construit : ::
143         
144                 extMesh = mc.MEDCouplingMappedExtrudedMesh(mesh3D, mesh2D, 0)
145         
146         On a alors la garantie que, dans ``extMesh``,  les cellules sont ordonnées par niveau Z croissant. 
147         Il suffit de récupérer le 2ème niveau (``MEDCouplingMappedExtrudedMesh.getMesh3DIds()``). ::
148         
149                 n_cells = mesh2D.getNumberOfCells()
150                 cellIdsSol3 = extMesh.getMesh3DIds()[n_cells:2*n_cells]
151
152 On vérifie alors que les 3 solutions sont les mêmes : ::
153
154         print cellIdsSol1.getValues()
155         print cellIdsSol2.getValues()
156         print cellIdsSol3.getValues()
157
158
159 Extraire une sous partie d'un maillage 3D
160 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
161
162 Utiliser les identifiants de cellules ``cellIdsSol2`` obtenus précédemment pour extraire une sous-partie de ``mesh3D``,
163 c'est-à-dire un maillage avec un sous-ensemble des cellules de ``mesh3D``. ::
164
165         mesh3DPart = mesh3D[cellIdsSol2] 
166         
167 .. note:: En C++ la méthode sous-jacente invoquée (et par ailleurs aussi disponible en Python) s'appelle    
168         ``mesh3DPart = mesh3D.buildPartOfMySelf(cellIdsSol2,True)``
169
170 .. note:: Le type géométrique ne rentre pas du tout en compte ici. L'instruction précédente prend les cellules
171         dans l'ordre où elles sont disponibles dans le maillage initial. 
172
173 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 
174 tableau new-2-old.
175
176 A ce point, ``mesh3DPart`` repose sur une copie du tableau de coordonnées de ``mesh3D``, c'est-à-dire  60 nodes. 
177 Seuls 30 sont effectivement utilisés.
178 Pour retirer les noeuds orphelins de ``mesh3DPart`` invoquer simplement ``MEDCouplingUMesh.zipCoords()``. ::
179
180         mesh3DPart.zipCoords()
181
182 Maintenant, ``mesh3DPart`` repose sur 30 nodes et possède 6 cellules. Pour être prêt aux I/O MED-fichier, il est 
183 alors important de voir si ``mesh3DPart`` est bien ordonné, c'est-à-dire si ses cellules sont bien rangées par type géométrique.
184 On commence par inspecter l'état actuel : ::
185
186         print mesh3DPart.advancedRepr()
187         
188 La fonction suivante fait le même travail : ::
189
190         print mesh3DPart.checkConsecutiveCellTypesAndOrder([mc.NORM_HEXA8, mc.NORM_POLYHED])
191
192 Ou bien : ::
193
194         print mesh3DPart.checkConsecutiveCellTypes()
195
196 On voit que ``mesh3DPart`` contient 6 cellules, quatre HEXA8 puis deux POLYHED. Les cellules sont bien 
197 groupées par type géométrique. Si ce n'était pas le cas, on aurait pu invoquer ``MEDCouplingUMesh.sortCellsInMEDFileFrmt()``.
198
199
200 Extraire des cellules alignées sur une ligne 3D
201 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
202
203 On souhaite extraire de ``mesh3D`` les 3 cellules dont les barycentres sont le long de la ligne portée par
204 ``v = [0.,0.,1.]`` et passant par ``pt = [250.,150.,0.]``.
205 Il y a deux solutions.
206
207 * les barycentres de ``mesh3D``  : même principe qu'au-dessus. ::
208
209         baryXY = bary[:,[0,1]]
210         baryXY -= [250.,150.]
211         magn = baryXY.magnitude()
212         cellIds2Sol1 = magn.findIdsInRange(0.,1e-12)
213         
214 * utiliser le maillage extrudé ``extMesh`` : partant de l'unique cellule dans ``mesh2D`` dont le centre est 
215   en ``[250.,150.,0.]``, la méthdode ``MEDCouplingMappedExtrudedMesh.getMesh3DIds()`` retourne les identifiants de 
216   cellules rangée par rangée. ::
217
218         bary2 = mesh2D.computeCellCenterOfMass()[:,[0,1]]
219         bary2 -= [250.,150.]
220         magn = bary2.magnitude()
221         ids = magn.findIdsInRange(0.,1e-12)
222         idStart = int(ids) # ids is assumed to contain only one value, if not an exception is thrown
223         ze_range = range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells())
224         cellIds2Sol2 = extMesh.getMesh3DIds()[ze_range]
225
226 Maintenant on construit cette sous partie de ``mesh3D`` en utilisant ``cellIds2Sol1`` ou ``cellIds2Sol2``: ::
227
228         mesh3DSlice2 = mesh3D[cellIds2Sol1]
229         mesh3DSlice2.zipCoords()
230
231 Duplication, translation et aggrégation de maillages
232 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
233
234 Cette partie de l'exercice est intéressante pour construire des maillages complexes, ou pour aggréger des parties 
235 de maillages venant de différents processeurs.
236
237 On cherche ici à dupliquer ``mesh3DSlice2``, le translater et l'aggréger avec l'original.
238
239 Effectuer une copie complète de ``mesh3DSlice2`` (aussi appelée *deep copy*) sous le nom ``mesh3DSlice2bis``. 
240 Sur cette copie effectuer une translation de ``v=[0.,1000.,0.]``.
241 Puis aggréger ``mesh3DSlice2`` avec sa copie translatée ``mesh3DSlice2bis``, en utilisant ``MEDCouplingUMesh.MergeUMeshes()``. ::
242
243         mesh3DSlice2bis = mesh3DSlice2.deepCopy()
244         mesh3DSlice2bis.translate([0.,1000.,0.])
245         mesh3DSlice2All = mc.MEDCouplingUMesh.MergeUMeshes([mesh3DSlice2,mesh3DSlice2bis])
246         mesh3DSlice2All.writeVTK("mesh3DSlice2All.vtu")
247
248 .. note:: Pour information pour merger deux (ou plus) maillages non structurés, il faut invoquer ``MEDCouplingUMesh.MergeUMeshes()``
249         puis ``MEDCouplingUMesh.mergeNodes()`` sur le résultat, et enfin ``MEDCouplingUMesh.zipConnectivityTraducer()``.
250
251 .. _exo-umesh-desc-connec:
252
253 Connectivité descendante
254 ~~~~~~~~~~~~~~~~~~~~~~~~
255
256 Le but ici est de présenter la notion de *connectivité descendante* (*descending connectivity*).
257
258 La connectivité descendante représente les éléments de dimension ``n-1`` 
259 constituant chacune des cellules de dimension ``n`` (avec donc ``n`` la dimension du maillage, *mesh-dim*). Par exemple, pour un
260 maillage de dimension 3 (les cellules sont des *volumes* 3D), cela donne l'ensemble des faces (des *surfaces* 2D) bordant
261 ces volumes.  
262
263 A titre d'exemple, on se propose dans notre cas de récupérer les faces *internes* du maillage ``mesh3D``.
264 Pour cela il est nécessaire de construire le maillage 
265 descendant de ``mesh3D`` (stocké dans ``mesh3DSurf``) c'est-à-dire 
266 le maillage de mesh-dimension 2 (soit ``mesh3D.getMeshDimension()-1``) constitué
267 des *faces* bordant chacune des cellules (ici des *volumes* 3D) de ``mesh3D``.
268 La méthode ``MEDCoupling.buildDescendingConnectivity()`` calcule ce maillage, et retourne en même temps des tableaux 
269 de correspondance. Ces tableaux font le lien entre les identifiants des cellules de ``mesh3D`` 
270 vers les identifiants de cellules de ``mesh3DSurf``, et vice-et-versa.
271
272 Une face de ``mesh3DSurf`` est dite interne, si et seulement si, elle est partagée par plus d'une cellule 3D de ``mesh3D``. 
273 Les 3ème et 4ème paramètres de sortie de la fonction donnent le lien 
274 entre une face et ses cellules *parentes* (i.e. le ou les volumes qu'elle délimite). 
275 Ce lien est exprimé au format *indirect index* vu dans le premier exercice :ref:`indirect-index-exo`. ::
276
277         mesh3DSurf, desc, descIndx, revDesc, revDescIndx = mesh3D.buildDescendingConnectivity()
278         numberOf3DCellSharing = revDescIndx.deltaShiftIndex()
279         cellIds = numberOf3DCellSharing.findIdsNotEqual(1)
280         mesh3DSurfInside = mesh3DSurf[cellIds]
281         mesh3DSurfInside.writeVTK("mesh3DSurfInside.vtu")
282         
283 Ce genre de manipulation est très utile pour accéder au voisinage d'une ou plusieurs cellules d'un maillage non-structuré. 
284  
285 .. image:: images/mesh3DSurfInside.jpg
286
287 Solution
288 ~~~~~~~~
289
290 :ref:`python_testMEDCouplingumesh1_solution`