2 Spliter et fusionner un fichier MED grâce à l'API avancée de MEDLoader
3 ----------------------------------------------------------------------
8 Cet exercise présente un cas complet et avancé d'utilisation de l'API avancée de MEDLoader.
9 Le but est de créer un maillage multi-type à partir de rien avec 2 champs :
11 * un champ aux cellules "CellField"
12 * un champ aux noeuds "NodeField"
14 Nous allons ensuite couper ces champs en deux parties (dans le but d'un traitement en parallèle par un code par exemple)
15 et aussi montrer comment re-fusionner deux champs à partir de morceaux disjoints.
17 Début de l'implémentation
18 ~~~~~~~~~~~~~~~~~~~~~~~~~
20 Créer un unstructured mesh ``m0`` issu d'un maillage structuré (meshDim=2, spaceDim=2) de 30*30.
21 Chacune des cellules paires du maillage sera *simplexisée* (i.e. coupée en triangle - méthode ``MEDCouplingUMesh.simplexize(0)``) ::
23 import MEDLoader as ml
24 from MEDLoader import MEDLoader
26 m0 = ml.MEDCouplingCMesh()
27 arr = ml.DataArrayDouble(31,1) ; arr.iota(0.)
29 m0 = m0.buildUnstructured()
30 m00 = m0[::2] # Extract even cells
33 m0 = ml.MEDCouplingUMesh.MergeUMeshes([m00,m01])
34 m0.getCoords()[:] *= 1/15.
37 .. note:: Le ``setName()`` sur "m0" est obligatoire. Ne pas oublier que dans le contexte MED fichier
38 le nommage correct des maillages est fondamental.
40 Créer les champs ``cellField`` et ``nodeField`` au pas de temps identifié à (5,6) et au pas de temps 5.6 s. ::
43 cellField = ml.MEDCouplingFieldDouble(ml.ON_CELLS, ml.ONE_TIME)
44 cellField.setTime(5.6,5,6)
46 cellField.setName("CellField")
47 cellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
48 cellField.getArray().setInfoOnComponent(0,"powercell [W]")
50 nodeField = ml.MEDCouplingFieldDouble(ml.ON_NODES,ml.ONE_TIME)
51 nodeField.setTime(5.6,5,6)
53 nodeField.setName("NodeField")
54 nodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
55 nodeField.getArray().setInfoOnComponent(0,"powernode [W]")
57 On obtient par exemple pour "CellField" ceci :
59 .. image:: images/SplitAndMergeCell1.jpg
62 Partitionnement de maillage
63 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
65 Couper ``m0`` en deux parties distinctes. Les deux parties seront nommées ``proc0`` et ``proc1``.
66 ``proc0`` sera la partie dans la boîte englobante (``MEDCouplingUMesh.getCellsInBoundingBox()``) ``[(0.,0.4),(0.,0.4)]``
67 à 1e-10 près. ``proc1`` sera le complémentaire (``DataArrayInt.buildComplement()``). ::
69 proc0 = m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
70 proc1 = proc0.buildComplement(m0.getNumberOfCells())
72 .. image:: images/SplitAndMerge2.jpg
74 Ecriture dans 2 fichiers MED séparés
75 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77 En partant du partitionnement ``proc0`` et ``proc1`` créer 2 fichiers MED appelés "proc0.med" et "proc1.med" : ::
79 nodeField0 = nodeField[proc0] ; cellField0 = cellField[proc0] ; cellField0.setMesh(nodeField0.getMesh())
80 nodeField1 = nodeField[proc1] ; cellField1 = cellField[proc1] ; cellField1.setMesh(nodeField1.getMesh())
82 proc0_fname = "proc0.med"
83 MEDLoader.WriteField(proc0_fname, nodeField0, True)
84 MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc0_fname, cellField0)
86 proc1_fname = "proc1.med"
87 MEDLoader.WriteField(proc1_fname,nodeField1,True)
88 MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,cellField1)
90 Lecture et fusion des 2 fichiers MED séparés (non optimal)
91 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
93 Partant de "proc0.med" et de "proc1.med" lire leur "CellField" respectif avec l'API basique,
94 agréger les deux et mettre le résultat dans ``cellField_read`` : ::
96 cellField0_read = MEDLoader.ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
97 cellField1_read = MEDLoader.ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
98 cellField_read = ml.MEDCouplingFieldDouble.MergeFields([cellField0_read,cellField1_read])
100 .. note:: On peut avoir l'impression que l'information Cell (méthode ``ReadFieldCell``) est répétée de manière abusive
101 (effectivement le champ "CellField" a été créé aux cellules),
102 mais ne pas oublier que dans la norme MED fichier rien n'interdit qu'un champ repose sur des cellules mais
103 aussi simultanément sur des noeuds, ou des points de Gauss ...
105 Comparer ``cellField_read`` et ``cellField0``. Problème, à cause de la contrainte sur la numérotation MED fichier,
106 on a perdu la numérotation initiale. Ou plus exactement il n'y a pas
107 de moyen standard de retrouver la numérotation originale. Donc un ``MEDCouplingFieldDouble.isEqual()``
108 n'est pas suffisant. Utilisons un ``MEDCouplingFieldDouble.substractInPlaceDM()``
109 qui opère pour nous une renumérotation suivant une politique particulière (*policy*, voir doc html).
110 Pour ce faire une copie profonde (*deep copy*) de ``cellField`` vers ``cellFieldCpy`` et opérer sur cette copie
111 un ``substractInPlaceDM`` (DM pour "Different Meshes", contrairement à ``substract`` qui ne marche que
112 s'ils partagent le même maillage): ::
114 cellFieldCpy = cellField.deepCpy()
115 cellFieldCpy.substractInPlaceDM(cellField_read,10,1e-12)
116 cellFieldCpy.getArray().abs()
117 print cellFieldCpy.getArray().isUniform(0.,1e-12)
119 Opérons le même travail sur "NodeField" que celui réalisé plus haut sur "CellField".
120 La différence ici c'est qu'il va y avoir duplication de l'information à la frontière, car les noeuds limites sont partagés
123 nodeField0_read = MEDLoader.ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
124 nodeField1_read = MEDLoader.ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
125 nodeField_read = ml.MEDCouplingFieldDouble.MergeFields([nodeField0_read, nodeField1_read])
127 .. note:: Dans cette partie, on a donc relu le maillage une deuxième fois ce qui peut être pénalisant ...
129 Invoquer ``MEDCouplingUMesh.mergeNodes()`` sur ``nodeField_read`` pour lui retirer les noeuds dupliqués.
130 Faire une deep copy appelée ``nodeFieldCpy`` de ``nodeField``
131 et supprimer encore les doublons : ::
133 nodeField_read.mergeNodes(1e-10)
134 nodeFieldCpy = nodeField.deepCpy()
135 nodeFieldCpy.mergeNodes(1e-10)
137 .. note:: A noter que ``mergeNodes()`` possède deux paramètres de précisions (*epsilons*), le premier,
138 classique, sur la distance absolue entre les noeuds, et l'autre sur la tolérance acceptée sur les valeurs du champ.
139 Si la valeur du champ de deux noeuds à fusionner dépasse ce deuxième epsilon, une exception est levée.
141 Comparer ``nodeFieldCpy`` et ``nodeField_read`` toujours en utilisant ``MEDCouplingFieldDouble.substractInPlaceDM()`` : ::
143 nodeFieldCpy.substractInPlaceDM(nodeField_read,10,1e-12)
144 print nodeFieldCpy.getArray().isUniform(0.,1e-12)
147 Lecture et merge des 2 fichiers MED séparés (moins facile, mais plus optimal)
148 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
150 Il s'agit ici de faire une méthode plus systématique et potentiellement plus générale de fusion de fichiers.
151 Pour de gros fichiers cette approche est à préférer.
152 Outre la performance, cette approche a l'avantage de pouvoir rajouter des infos.
154 Avec l'API avancée lire les maillages des deux fichiers "proc0.med" et "proc1.med" et agréger le résultat
155 dans une instance ``mergeMLMesh`` de ``MEDFileUMesh``.
156 Traiter tous les niveaux de dimension (même si ici il n'y en a qu'un seul) en utilisant la méthode ``MEDFileUMesh.getNonEmptyLevels()``
157 sur l'instance venant de "proc0.med".
159 La solution donnée ci-dessous est la plus générique possible, car elle traite aussi les différents pas de temps et les
160 différents types géométriques : ::
162 fileNames = ["proc0.med","proc1.med"]
163 msML = [ml.MEDFileMesh.New(fname) for fname in fileNames]
164 fsML = [ml.MEDFileFields.New(fname) for fname in fileNames]
165 mergeMLMesh = ml.MEDFileUMesh()
166 mergeMLFields = ml.MEDFileFields()
167 for lev in msML[0].getNonEmptyLevels():
168 o2nML = len(msML[0].getNonEmptyLevels())*[None]
169 cs = [mML.getCoords() for mML in msML]
170 mergeMLMesh.setCoords(ml.DataArrayDouble.Aggregate(cs))
171 ms = [mML.getMeshAtLevel(lev) for mML in msML]
172 m = ml.MEDCouplingUMesh.MergeUMeshes(ms) ; m.setCoords(mergeMLMesh.getCoords())
173 o2nML[lev] = m.sortCellsInMEDFileFrmt()
174 mergeMLMesh.setMeshAtLevel(lev,m)
177 for fieldName in fsML[0].getFieldsNames():
178 fmts = [fML[fieldName] for fML in fsML]
179 mergeField = ml.MEDFileFieldMultiTS()
180 for dt,it,tim in fmts[0].getTimeSteps():
181 fts = [fmt[dt,it] for fmt in fmts]
182 arrs = len(fts)*[None]
183 for typp in fts[0].getTypesOfFieldAvailable():
185 if typp == ml.ON_CELLS:
187 for geoTyp,smth in ft.getFieldSplitedByType():
188 if geoTyp != ml.NORM_ERROR:
189 smth1 = filter(lambda x:x[0] == ml.ON_CELLS,smth)
190 arr2s = [ft.getUndergroundDataArray()[elt[1][0]:elt[1][1]] for elt in smth1]
191 arr1s.append(ml.DataArrayDouble.Aggregate(arr2s))
198 smth = filter(lambda x:x[0] == ml.NORM_ERROR,ft.getFieldSplitedByType())
199 arr2 = ml.DataArrayDouble.Aggregate([ft.getUndergroundDataArray()[elt[1][0][1][0]:elt[1][0][1][1]] for elt in smth])
203 arr = ml.DataArrayDouble.Aggregate(arr1s)
204 if typp == ml.ON_CELLS:
205 arr.renumberInPlace(o2nML[lev])
206 mcf = ml.MEDCouplingFieldDouble(typp,ml.ONE_TIME) ; mcf.setName(fieldName) ; mcf.setTime(tim,dt,it) ; mcf.setArray(arr)
207 mcf.setMesh(mergeMLMesh.getMeshAtLevel(lev)) ; mcf.checkCoherency()
208 mergeField.appendFieldNoProfileSBT(mcf)
211 mergeMLFields.pushField(mergeField)
213 mergeMLMesh.write("merge.med",2)
214 mergeMLFields.write("merge.med",0)
220 :ref:`python_testMEDLoaderSplitAndMerge1_solution`