doc/tutorial/atestMEDLoaderSplitAndMerge1.rst

   1
   2 .. _python_testMEDLoaderSplitAndMerge1_solution:
   3
   4 Splitting and Merging a MED file using MEDLoader
   5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   6
   7 ::
   8
   9         import MEDLoader as ml
  10         from MEDLoader import MEDLoader
  11
  12         m0 = ml.MEDCouplingCMesh()
  13         arr = ml.DataArrayDouble(31,1) ; arr.iota(0.)
  14         m0.setCoords(arr,arr)
  15         m0 = m0.buildUnstructured()
  16         m00 = m0[::2]      # Extract even cells
  17         m00.simplexize(0)
  18         m01 = m0[1::2]
  19         m0 = ml.MEDCouplingUMesh.MergeUMeshes([m00,m01])
  20         m0.getCoords()[:] *= 1/15.
  21         m0.setName("mesh")
  22         # Cell field
  23         cellField = ml.MEDCouplingFieldDouble(ml.ON_CELLS, ml.ONE_TIME)
  24         cellField.setTime(5.6,5,6)
  25         cellField.setMesh(m0)
  26         cellField.setName("CellField")
  27         cellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
  28         cellField.getArray().setInfoOnComponent(0,"powercell [W]")
  29         # Node field
  30         nodeField = ml.MEDCouplingFieldDouble(ml.ON_NODES,ml.ONE_TIME)
  31         nodeField.setTime(5.6,5,6)
  32         nodeField.setMesh(m0)
  33         nodeField.setName("NodeField")
  34         nodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
  35         nodeField.getArray().setInfoOnComponent(0,"powernode [W]")
  36         # Splitting
  37         proc0 = m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
  38         proc1 = proc0.buildComplement(m0.getNumberOfCells())
  39         #
  40         nodeField0 = nodeField[proc0] ; cellField0 = cellField[proc0] ; cellField0.setMesh(nodeField0.getMesh())
  41         nodeField1 = nodeField[proc1] ; cellField1 = cellField[proc1] ; cellField1.setMesh(nodeField1.getMesh())
  42
  43         proc0_fname = "proc0.med"
  44         MEDLoader.WriteField(proc0_fname, nodeField0, True)
  45         MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc0_fname, cellField0)
  46
  47         proc1_fname = "proc1.med"
  48         MEDLoader.WriteField(proc1_fname,nodeField1,True)
  49         MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,cellField1)
  50         #
  51         # Merging - Sub-optimal method
  52         #
  53         cellField0_read = MEDLoader.ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
  54         cellField1_read = MEDLoader.ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
  55         cellField_read = ml.MEDCouplingFieldDouble.MergeFields([cellField0_read,cellField1_read])
  56         cellFieldCpy = cellField.deepCpy()
  57         cellFieldCpy.substractInPlaceDM(cellField_read,10,1e-12)
  58         cellFieldCpy.getArray().abs()
  59         print cellFieldCpy.getArray().isUniform(0.,1e-12)
  60         #
  61         nodeField0_read = MEDLoader.ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
  62         nodeField1_read = MEDLoader.ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
  63         nodeField_read = ml.MEDCouplingFieldDouble.MergeFields([nodeField0_read, nodeField1_read])
  64         nodeField_read.mergeNodes(1e-10)
  65         nodeFieldCpy = nodeField.deepCpy()
  66         nodeFieldCpy.mergeNodes(1e-10)
  67         nodeFieldCpy.substractInPlaceDM(nodeField_read,10,1e-12)
  68         print nodeFieldCpy.getArray().isUniform(0.,1e-12)
  69         #
  70         # Merging - Optimal method
  71         #
  72         fileNames = ["proc0.med","proc1.med"]
  73         msML = [ml.MEDFileMesh.New(fname) for fname in fileNames]
  74         fsML = [ml.MEDFileFields.New(fname) for fname in fileNames]
  75         mergeMLMesh = ml.MEDFileUMesh()
  76         mergeMLFields = ml.MEDFileFields()
  77         for lev in msML[0].getNonEmptyLevels():
  78                 o2nML = len(msML[0].getNonEmptyLevels())*[None]
  79                 cs = [mML.getCoords() for mML in msML]
  80                 mergeMLMesh.setCoords(ml.DataArrayDouble.Aggregate(cs))
  81                 ms = [mML.getMeshAtLevel(lev) for mML in msML]
  82                 m = ml.MEDCouplingUMesh.MergeUMeshes(ms) ; m.setCoords(mergeMLMesh.getCoords())
  83                 o2nML[lev] = m.sortCellsInMEDFileFrmt()
  84                 mergeMLMesh.setMeshAtLevel(lev,m)
  85                 pass
  86
  87         for fieldName in fsML[0].getFieldsNames():
  88                 fmts = [fML[fieldName] for fML in fsML]
  89                 mergeField = ml.MEDFileFieldMultiTS()
  90                 for dt,it,tim in fmts[0].getTimeSteps():
  91                         fts = [fmt[dt,it] for fmt in fmts]
  92                         arrs = len(fts)*[None]
  93                         for typp in fts[0].getTypesOfFieldAvailable():
  94                                 arr1s = []
  95                                 if typp == ml.ON_CELLS:
  96                                         for ft in fts:
  97                                                 for geoTyp,smth in ft.getFieldSplitedByType():
  98                                                         if geoTyp != ml.NORM_ERROR:
  99                                                                 smth1 = filter(lambda x:x[0] == ml.ON_CELLS,smth)
 100                                                                 arr2s = [ft.getUndergroundDataArray()[elt[1][0]:elt[1][1]] for elt in smth1]
 101                                                                 arr1s.append(ml.DataArrayDouble.Aggregate(arr2s))
 102                                                                 pass
 103                                                         pass
 104                                                 pass
 105                                         pass
 106                                 else:
 107                                         for ft in fts:
 108                                                 smth = filter(lambda x:x[0] == ml.NORM_ERROR,ft.getFieldSplitedByType())
 109                                                 arr2 = ml.DataArrayDouble.Aggregate([ft.getUndergroundDataArray()[elt[1][0][1][0]:elt[1][0][1][1]] for elt in smth])
 110                                                 arr1s.append(arr2)
 111                                                 pass
 112                                         pass
 113                                 arr = ml.DataArrayDouble.Aggregate(arr1s)
 114                                 if typp == ml.ON_CELLS:
 115                                      arr.renumberInPlace(o2nML[lev])
 116                                 mcf = ml.MEDCouplingFieldDouble(typp,ml.ONE_TIME) ; mcf.setName(fieldName) ; mcf.setTime(tim,dt,it) ; mcf.setArray(arr)
 117                                 mcf.setMesh(mergeMLMesh.getMeshAtLevel(lev)) ; mcf.checkCoherency()
 118                                 mergeField.appendFieldNoProfileSBT(mcf)
 119                                 pass
 120                         pass
 121                 mergeMLFields.pushField(mergeField)
 122                 pass
 123         mergeMLMesh.write("merge.med",2)
 124         mergeMLFields.write("merge.med",0)