doc/tutorial/medloader_SplitAndMerge1_fr.rst

   1
   2 Spliter et fusionner un fichier MED grâce à l'API avancée de MEDLoader
   3 ----------------------------------------------------------------------
   4
   5 Objectif
   6 ~~~~~~~~
   7
   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 :
  10
  11 * un champ aux cellules "CellField"
  12 * un champ aux noeuds "NodeField"
  13
  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.
  16
  17 Début de l'implémentation
  18 ~~~~~~~~~~~~~~~~~~~~~~~~~
  19
  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)``) ::
  22
  23         import MEDLoader as ml
  24         from MEDLoader import MEDLoader
  25
  26         m0 = ml.MEDCouplingCMesh()
  27         arr = ml.DataArrayDouble(31,1) ; arr.iota(0.)
  28         m0.setCoords(arr,arr)
  29         m0 = m0.buildUnstructured()
  30         m00 = m0[::2]      # Extract even cells
  31         m00.simplexize(0)
  32         m01 = m0[1::2]
  33         m0 = ml.MEDCouplingUMesh.MergeUMeshes([m00,m01])
  34         m0.getCoords()[:] *= 1/15.
  35         m0.setName("mesh")
  36
  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.
  39
  40 Créer les champs ``cellField`` et ``nodeField`` au pas de temps identifié à (5,6) et au pas de temps 5.6 s. ::
  41
  42         # Cell field
  43         cellField = ml.MEDCouplingFieldDouble(ml.ON_CELLS, ml.ONE_TIME)
  44         cellField.setTime(5.6,5,6)
  45         cellField.setMesh(m0)
  46         cellField.setName("CellField")
  47         cellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
  48         cellField.getArray().setInfoOnComponent(0,"powercell [W]")
  49         # Node field
  50         nodeField = ml.MEDCouplingFieldDouble(ml.ON_NODES,ml.ONE_TIME)
  51         nodeField.setTime(5.6,5,6)
  52         nodeField.setMesh(m0)
  53         nodeField.setName("NodeField")
  54         nodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
  55         nodeField.getArray().setInfoOnComponent(0,"powernode [W]")
  56
  57 On obtient par exemple pour "CellField" ceci :
  58
  59 .. image:: images/SplitAndMergeCell1.jpg
  60
  61
  62 Partitionnement de maillage
  63 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
  64
  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()``). ::
  68
  69         proc0 = m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
  70         proc1 = proc0.buildComplement(m0.getNumberOfCells())
  71
  72 .. image:: images/SplitAndMerge2.jpg
  73
  74 Ecriture dans 2 fichiers MED séparés
  75 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  76
  77 En partant du partitionnement ``proc0`` et ``proc1`` créer 2 fichiers MED appelés "proc0.med" et "proc1.med" : ::
  78
  79         nodeField0 = nodeField[proc0] ; cellField0 = cellField[proc0] ; cellField0.setMesh(nodeField0.getMesh())
  80         nodeField1 = nodeField[proc1] ; cellField1 = cellField[proc1] ; cellField1.setMesh(nodeField1.getMesh())
  81
  82         proc0_fname = "proc0.med"
  83         MEDLoader.WriteField(proc0_fname, nodeField0, True)
  84         MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc0_fname, cellField0)
  85
  86         proc1_fname = "proc1.med"
  87         MEDLoader.WriteField(proc1_fname,nodeField1,True)
  88         MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,cellField1)
  89
  90 Lecture et fusion des 2 fichiers MED séparés (non optimal)
  91 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  92
  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`` : ::
  95
  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])
  99
 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 ...
 104
 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): ::
 113
 114         cellFieldCpy = cellField.deepCpy()
 115         cellFieldCpy.substractInPlaceDM(cellField_read,10,1e-12)
 116         cellFieldCpy.getArray().abs()
 117         print cellFieldCpy.getArray().isUniform(0.,1e-12)
 118
 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
 121 des deux côtés : ::
 122
 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])
 126
 127 .. note:: Dans cette partie, on a donc relu le maillage une deuxième fois ce qui peut être pénalisant ...
 128
 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 : ::
 132
 133         nodeField_read.mergeNodes(1e-10)
 134         nodeFieldCpy = nodeField.deepCpy()
 135         nodeFieldCpy.mergeNodes(1e-10)
 136
 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.
 140
 141 Comparer ``nodeFieldCpy`` et ``nodeField_read`` toujours en utilisant ``MEDCouplingFieldDouble.substractInPlaceDM()`` : ::
 142
 143         nodeFieldCpy.substractInPlaceDM(nodeField_read,10,1e-12)
 144         print nodeFieldCpy.getArray().isUniform(0.,1e-12)
 145
 146
 147 Lecture et merge des 2 fichiers MED séparés (moins facile, mais plus optimal)
 148 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 149
 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.
 153
 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".
 158
 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 : ::
 161
 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)
 175                 pass
 176
 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():
 184                                 arr1s = []
 185                                 if typp == ml.ON_CELLS:
 186                                         for ft in fts:
 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))
 192                                                                 pass
 193                                                         pass
 194                                                 pass
 195                                         pass
 196                                 else:
 197                                         for ft in fts:
 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])
 200                                                 arr1s.append(arr2)
 201                                                 pass
 202                                         pass
 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)
 209                                 pass
 210                         pass
 211                 mergeMLFields.pushField(mergeField)
 212                 pass
 213         mergeMLMesh.write("merge.med",2)
 214         mergeMLFields.write("merge.med",0)
 215
 216
 217 Solution
 218 ~~~~~~~~
 219
 220 :ref:`python_testMEDLoaderSplitAndMerge1_solution`