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 medcoupling as mc
  24
  25         m0 = mc.MEDCouplingCMesh()
  26         arr = mc.DataArrayDouble(31,1) ; arr.iota(0.)
  27         m0.setCoords(arr,arr)
  28         m0 = m0.buildUnstructured()
  29         m00 = m0[::2]                # Extract even cells
  30         m00.simplexize(0)
  31         m01 = m0[1::2]
  32         m0 = mc.MEDCouplingUMesh.MergeUMeshes([m00,m01])
  33         m0.getCoords()[:] *= 1/15.   # Illustrate how to quickly rescale a mesh
  34         m0.setName("mesh")
  35
  36 .. note:: Le ``setName()`` sur "m0" est obligatoire. Ne pas oublier que dans le contexte MED fichier
  37         le nommage correct des maillages est fondamental.
  38
  39 Créer les champs ``cellField`` et ``nodeField`` au pas de temps identifié à (5,6) et au pas de temps 5.6 s. ::
  40
  41         # Cell field
  42         cellField = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
  43         cellField.setTime(5.6,5,6)
  44         cellField.setMesh(m0)
  45         cellField.setName("CellField")
  46         cellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
  47         cellField.getArray().setInfoOnComponent(0,"powercell [W]")
  48         # Node field
  49         nodeField = mc.MEDCouplingFieldDouble(mc.ON_NODES,mc.ONE_TIME)
  50         nodeField.setTime(5.6,5,6)
  51         nodeField.setMesh(m0)
  52         nodeField.setName("NodeField")
  53         nodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
  54         nodeField.getArray().setInfoOnComponent(0,"powernode [W]")
  55
  56 On obtient par exemple pour "CellField" ceci :
  57
  58 .. image:: images/SplitAndMergeCell1.jpg
  59
  60
  61 Partitionnement de maillage
  62 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
  63
  64 Couper ``m0`` en deux parties distinctes. Les deux parties seront nommées ``proc0`` et ``proc1``.
  65 ``proc0`` sera la partie dans la boîte englobante (``MEDCouplingUMesh.getCellsInBoundingBox()``) ``[(0.,0.4),(0.,0.4)]``
  66 à 1e-10 près. ``proc1`` sera le complémentaire (``DataArrayInt.buildComplement()``). ::
  67
  68         proc0 = m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
  69         proc1 = proc0.buildComplement(m0.getNumberOfCells())
  70
  71 .. image:: images/SplitAndMerge2.jpg
  72
  73 Ecriture dans 2 fichiers MED séparés
  74 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  75
  76 En partant du partitionnement ``proc0`` et ``proc1`` créer 2 fichiers MED appelés "proc0.med" et "proc1.med" : ::
  77
  78         nodeField0 = nodeField[proc0] ; cellField0 = cellField[proc0] ; cellField0.setMesh(nodeField0.getMesh())
  79         nodeField1 = nodeField[proc1] ; cellField1 = cellField[proc1] ; cellField1.setMesh(nodeField1.getMesh())
  80
  81         proc0_fname = "proc0.med"
  82         mc.WriteField(proc0_fname, nodeField0, True)
  83         mc.WriteFieldUsingAlreadyWrittenMesh(proc0_fname, cellField0)
  84
  85         proc1_fname = "proc1.med"
  86         mc.WriteField(proc1_fname,nodeField1,True)
  87         mc.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,cellField1)
  88
  89 Lecture et fusion des 2 fichiers MED séparés (non optimal)
  90 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  91
  92 Partant de "proc0.med" et de "proc1.med" lire leur "CellField" respectif avec l'API basique,
  93 agréger les deux et mettre le résultat dans ``cellField_read`` : ::
  94
  95         cellField0_read = mc.ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
  96         cellField1_read = mc.ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
  97         cellField_read = mc.MEDCouplingFieldDouble.MergeFields([cellField0_read,cellField1_read])
  98
  99 .. note:: On peut avoir l'impression que l'information Cell (méthode ``ReadFieldCell``) est répétée de manière abusive
 100         (effectivement le champ "CellField" a été créé aux cellules),
 101         mais ne pas oublier que dans la norme MED fichier rien n'interdit qu'un champ repose sur des cellules mais
 102         aussi simultanément sur des noeuds, ou des points de Gauss ...
 103
 104 Comparer ``cellField_read`` et ``cellField0``. Problème, à cause de la contrainte sur la numérotation MED fichier,
 105 on a perdu la numérotation initiale. Ou plus exactement il n'y a pas
 106 de moyen standard de retrouver la numérotation originale. Donc un ``MEDCouplingFieldDouble.isEqual()``
 107 n'est pas suffisant. Utilisons un ``MEDCouplingFieldDouble.substractInPlaceDM()``
 108 qui opère pour nous une renumérotation suivant une politique particulière (*policy*, voir doc html).
 109 Pour ce faire une copie profonde (*deep copy*) de ``cellField`` vers ``cellFieldCpy`` et opérer sur cette copie
 110 un ``substractInPlaceDM`` (DM pour "Different Meshes", contrairement à ``substract`` qui ne marche que
 111 s'ils partagent le même maillage): ::
 112
 113         cellFieldCpy = cellField.deepCopy()
 114         cellFieldCpy.substractInPlaceDM(cellField_read,10,1e-12)
 115         cellFieldCpy.getArray().abs()
 116         print(cellFieldCpy.getArray().isUniform(0.,1e-12))
 117
 118 Opérons le même travail sur "NodeField" que celui réalisé plus haut sur "CellField".
 119 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
 120 des deux côtés : ::
 121
 122         nodeField0_read = mc.ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
 123         nodeField1_read = mc.ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
 124         nodeField_read = mc.MEDCouplingFieldDouble.MergeFields([nodeField0_read, nodeField1_read])
 125
 126 .. note:: Dans cette partie, on a donc relu le maillage une deuxième fois ce qui peut être pénalisant ...
 127
 128 Invoquer ``MEDCouplingUMesh.mergeNodes()`` sur ``nodeField_read`` pour lui retirer les noeuds dupliqués.
 129 Faire une deep copy appelée ``nodeFieldCpy`` de ``nodeField``
 130 et supprimer encore les doublons : ::
 131
 132         nodeField_read.mergeNodes(1e-10)
 133         nodeFieldCpy = nodeField.deepCopy()
 134         nodeFieldCpy.mergeNodes(1e-10)
 135
 136 .. note:: A noter que ``mergeNodes()`` possède deux paramètres de précisions (*epsilons*), le premier,
 137         classique, sur la distance absolue entre les noeuds, et l'autre sur la tolérance acceptée sur les valeurs du champ.
 138         Si la valeur du champ de deux noeuds à fusionner dépasse ce deuxième epsilon, une exception est levée.
 139
 140 Comparer ``nodeFieldCpy`` et ``nodeField_read`` toujours en utilisant ``MEDCouplingFieldDouble.substractInPlaceDM()`` : ::
 141
 142         nodeFieldCpy.substractInPlaceDM(nodeField_read,10,1e-12)
 143         print(nodeFieldCpy.getArray().isUniform(0.,1e-12))
 144
 145
 146 Lecture et merge des 2 fichiers MED séparés (moins facile, mais plus optimal)
 147 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 148
 149 Il s'agit ici de faire une méthode plus systématique et potentiellement plus générale de fusion de fichiers.
 150 Pour de gros fichiers cette approche est à préférer.
 151 Outre la performance, cette approche a l'avantage de pouvoir rajouter des infos.
 152
 153 Avec l'API avancée lire les maillages des deux fichiers "proc0.med" et "proc1.med" et agréger le résultat
 154 dans une instance ``mergeMLMesh`` de ``MEDFileUMesh``.
 155 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()``
 156 sur l'instance venant de "proc0.med".
 157
 158 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
 159 différents types géométriques : ::
 160
 161         fileNames = ["proc0.med","proc1.med"]
 162         msML = [mc.MEDFileMesh.New(fname) for fname in fileNames]
 163         fsML = [mc.MEDFileFields.New(fname) for fname in fileNames]
 164         mergeMLMesh = mc.MEDFileUMesh()
 165         mergeMLFields = mc.MEDFileFields()
 166         for lev in msML[0].getNonEmptyLevels():
 167                 o2nML = len(msML[0].getNonEmptyLevels())*[None]
 168                 cs = [mML.getCoords() for mML in msML]
 169                 mergeMLMesh.setCoords(mc.DataArrayDouble.Aggregate(cs))
 170                 ms = [mML.getMeshAtLevel(lev) for mML in msML]
 171                 m = mc.MEDCouplingUMesh.MergeUMeshes(ms) ; m.setCoords(mergeMLMesh.getCoords())
 172                 o2nML[lev] = m.sortCellsInMEDFileFrmt()
 173                 mergeMLMesh.setMeshAtLevel(lev,m)
 174                 pass
 175
 176         for fieldName in fsML[0].getFieldsNames():
 177                 fmts = [fML[fieldName] for fML in fsML]
 178                 mergeField = mc.MEDFileFieldMultiTS()
 179                 for dt,it,tim in fmts[0].getTimeSteps():
 180                         fts = [fmt[dt,it] for fmt in fmts]
 181                         arrs = len(fts)*[None]
 182                         for typp in fts[0].getTypesOfFieldAvailable():
 183                                 arr1s = []
 184                                 if typp == mc.ON_CELLS:
 185                                         for ft in fts:
 186                                                 for geoTyp,smth in ft.getFieldSplitedByType():
 187                                                         if geoTyp != mc.NORM_ERROR:
 188                                                                 smth1 = filter(lambda x:x[0] == mc.ON_CELLS,smth)
 189                                                                 arr2s = [ft.getUndergroundDataArray()[elt[1][0]:elt[1][1]] for elt in smth1]
 190                                                                 arr1s.append(mc.DataArrayDouble.Aggregate(arr2s))
 191                                                                 pass
 192                                                         pass
 193                                                 pass
 194                                         pass
 195                                 else:
 196                                         for ft in fts:
 197                                                 smth = filter(lambda x:x[0] == mc.NORM_ERROR,ft.getFieldSplitedByType())
 198                                                 arr2 = mc.DataArrayDouble.Aggregate([ft.getUndergroundDataArray()[elt[1][0][1][0]:elt[1][0][1][1]] for elt in smth])
 199                                                 arr1s.append(arr2)
 200                                                 pass
 201                                         pass
 202                                 arr = mc.DataArrayDouble.Aggregate(arr1s)
 203                                 if typp == mc.ON_CELLS:
 204                                      arr.renumberInPlace(o2nML[lev])
 205                                 mcf = mc.MEDCouplingFieldDouble(typp,mc.ONE_TIME) ; mcf.setName(fieldName) ; mcf.setTime(tim,dt,it) ; mcf.setArray(arr)
 206                                 mcf.setMesh(mergeMLMesh.getMeshAtLevel(lev)) ; mcf.checkConsistencyLight()
 207                                 mergeField.appendFieldNoProfileSBT(mcf)
 208                                 pass
 209                         pass
 210                 mergeMLFields.pushField(mergeField)
 211                 pass
 212         mergeMLMesh.write("merge.med",2)
 213         mergeMLFields.write("merge.med",0)
 214
 215
 216 Solution
 217 ~~~~~~~~
 218
 219 :ref:`python_testMEDLoaderSplitAndMerge1_solution`