doc/tutorial/medloader_SplitAndMerge1_en.rst

   1
   2 Splitting and merging a MED file using MEDLoader's advanced API
   3 ---------------------------------------------------------------
   4
   5 Objective
   6 ~~~~~~~~~
   7
   8 The aim of this exercise is to create a mesh with mixed geometrical type from scratch, and to associate two fields to it:
   9
  10 * a field on cells "CellField"
  11 * a field on nodes "NodeField"
  12
  13
  14 Implementation start
  15 ~~~~~~~~~~~~~~~~~~~~
  16
  17 Create an unstructured mesh "m0" built from a 30x30 structured mesh (meshDim=2, spaceDim=2).
  18 Each of the even cell of "m0" is "simplexized" (cut in triangles - method MEDCouplingUMesh.simplexize(0)) ::
  19
  20         from MEDLoader import *
  21         m0=MEDCouplingCMesh()
  22         arr=DataArrayDouble(31,1) ; arr.iota(0.)
  23         m0.setCoords(arr,arr)
  24         m0=m0.buildUnstructured()
  25         m00=m0[::2] ; m00.simplexize(0) ; m01=m0[1::2]
  26         m0=MEDCouplingUMesh.MergeUMeshes([m00,m01])
  27         m0.getCoords()[:]*=1/15.
  28         m0.setName("mesh")
  29
  30 .. note:: The call to setName() on "m0" is mandatory. Don't forget that the correct naming of the meshes is paramount in the MED file context.
  31
  32 Create the fields "CellField" and "NodeField" at the time-stamp (5,6) corresponding to 5.6s.
  33 ::
  34
  35         CellField=MEDCouplingFieldDouble(ON_CELLS,ONE_TIME) ; CellField.setTime(5.6,5,6) ; CellField.setMesh(m0)
  36         CellField.setName("CellField")
  37         CellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))") ; CellField.getArray().setInfoOnComponent(0,"powercell [W]")
  38         NodeField=MEDCouplingFieldDouble(ON_NODES,ONE_TIME) ; NodeField.setTime(5.6,5,6) ; NodeField.setMesh(m0)
  39         NodeField.setName("NodeField")
  40         NodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))") ; NodeField.getArray().setInfoOnComponent(0,"powernode [W]")
  41
  42 "CellField" looks like this:
  43
  44 .. image:: images/SplitAndMergeCell1.jpg
  45
  46
  47 Mesh partitionning
  48 ~~~~~~~~~~~~~~~~~~
  49
  50 Cut "m0" into two distinct parts called "proc0" and "proc1". "proc0" will be contained in the bounding box [(0.,0.4),(0.,0.4)] (with a precision of 1e-10). Use the method MEDCouplingUMesh.getCellsInBoundingBox(). "proc1" is simply the complementary part of "proc0" (method DataArrayInt.buildComplement()). ::
  51
  52      proc0=m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
  53      proc1=proc0.buildComplement(m0.getNumberOfCells())
  54
  55 .. image:: images/SplitAndMerge2.jpg
  56
  57 Writing into two different MED files
  58 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  59
  60 Starting with the partition above ("proc0" and "proc1") create two MED files called "proc0.med" et "proc1.med". ::
  61
  62      NodeField0=NodeField[proc0] ; CellField0=CellField[proc0] ; CellField0.setMesh(NodeField0.getMesh())
  63      NodeField1=NodeField[proc1] ; CellField1=CellField[proc1] ; CellField1.setMesh(NodeField1.getMesh())
  64
  65      proc0_fname="proc0.med"
  66      MEDLoader.WriteField(proc0_fname,NodeField0,True)
  67      MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc0_fname,CellField0)
  68
  69      proc1_fname="proc1.med"
  70      MEDLoader.WriteField(proc1_fname,NodeField1,True)
  71      MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,CellField1)
  72
  73 Reading and merging 2 MED files - Easy (but non optimal) version
  74 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  75
  76 In the two files "proc0.med" and "proc1.med" read the respective "CellField" with the basic API. Aggregate the two and store the result in "CellField_read". ::
  77
  78      CellField0_read=MEDLoader.ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
  79      CellField1_read=MEDLoader.ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
  80      CellField_read=MEDCouplingFieldDouble.MergeFields([CellField0_read,CellField1_read])
  81
  82 .. note:: It might seem to the reader that the cell type information is repeated uselessly, but don't forget that the MED file norm doesn't forbid a field to be defined simultaneously on nodes and on Gauss points for example ...
  83
  84 Compare "CellField_read" and "CellField0". Problem: because of the constraint on the MED file numbering, the initial numbering has been lost. Or more exactly there is no standard way to retrieve it. This means that a call to MEDCouplingFieldDouble.isEqual() won't succeed. Let's use the method MEDCouplingFieldDouble.substractInPlaceDM() which operates a renumbering based on a given policy (see HTML doc).
  85 To this end, create a deep copy of "CellField" into "CellFieldCpy" and invoke substractInPlaceDM() on it (DM stands for "Different Meshes", contrarily to substract() which only succeeds if the fields share the same mesh). ::
  86
  87      CellFieldCpy=CellField.deepCpy()
  88      CellFieldCpy.substractInPlaceDM(CellField_read,10,1e-12)
  89      CellFieldCpy.getArray().abs()
  90      print CellFieldCpy.getArray().isUniform(0.,1e-12)
  91
  92 Let's do the same on "NodeField". The main difference here is that redundant information is created at the boundary. ::
  93
  94      NodeField0_read=MEDLoader.ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
  95      NodeField1_read=MEDLoader.ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
  96      NodeField_read=MEDCouplingFieldDouble.MergeFields([NodeField0_read,NodeField1_read])
  97
  98 .. note:: The mesh is read a second time here, which can be damaging in terms of performance.
  99
 100 Invoke MEDCouplingUMesh.mergeNodes() on "NodeField_read" to remove duplicate nodes.
 101 Make a deep copy called  "NodeFieldCpy" from "NodeField" and call  MEDCouplingUMesh.mergeNodes(). ::
 102
 103      NodeField_read.mergeNodes(1e-10)
 104      NodeFieldCpy=NodeField.deepCpy()
 105      NodeFieldCpy.mergeNodes(1e-10)
 106
 107 .. note:: mergeNodes() takes two epsilons: the first classical one on the absolute distance between nodes, and the second expressing a tolerance on the values. If the field value of two nodes to be merged is bigger than this an exception is raised.
 108
 109 Compare "NodeFieldCpy" and "NodeField_read" still using MEDCouplingFieldDouble.substractInPlaceDM(). ::
 110
 111      NodeFieldCpy.substractInPlaceDM(NodeField_read,10,1e-12)
 112      print NodeFieldCpy.getArray().isUniform(0.,1e-12)
 113
 114
 115 Read/write of two separated MED files - More complex but more efficient version
 116 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 117
 118 We show a more systematic and more general method to merge files.
 119 This is the preferred route when dealing with big files .
 120 This method adds performance but also allows to add extra information.
 121
 122 Using the advanced API read the meshes of two files "proc0.med" and "proc1.med" and aggregate the result in an MEDFileUMesh instance "mergeMLMesh".
 123 Handle all the levels (even if there is only one in the present case) using the method
 124 MEDFileUMesh.getNonEmptyLevels() on the instance coming from "proc0.med". ::
 125
 126      fileNames=["proc0.med","proc1.med"]
 127      msML=[MEDFileMesh.New(fname) for fname in fileNames]
 128      fsML=[MEDFileFields.New(fname) for fname in fileNames]
 129      mergeMLMesh=MEDFileUMesh()
 130      mergeMLFields=MEDFileFields()
 131      for lev in msML[0].getNonEmptyLevels():
 132          o2nML=len(msML[0].getNonEmptyLevels())*[None]
 133          cs=[mML.getCoords() for mML in msML]
 134          mergeMLMesh.setCoords(DataArrayDouble.Aggregate(cs))
 135          ms=[mML.getMeshAtLevel(lev) for mML in msML]
 136          m=MEDCouplingUMesh.MergeUMeshes(ms) ; m.setCoords(mergeMLMesh.getCoords())
 137          o2nML[lev]=m.sortCellsInMEDFileFrmt()
 138          mergeMLMesh.setMeshAtLevel(lev,m)
 139          pass
 140
 141      for fieldName in fsML[0].getFieldsNames():
 142          fmts=[fML[fieldName] for fML in fsML]
 143          mergeField=MEDFileFieldMultiTS()
 144          for dt,it,tim in fmts[0].getTimeSteps():
 145              fts=[fmt[dt,it] for fmt in fmts]
 146              arrs=len(fts)*[None]
 147              for typp in fts[0].getTypesOfFieldAvailable():
 148                  arr1s=[]
 149                  if typp==ON_CELLS:
 150                      for ft in fts:
 151                          for geoTyp,smth in ft.getFieldSplitedByType():
 152                              if geoTyp!=NORM_ERROR:
 153                                  smth1=filter(lambda x:x[0]==ON_CELLS,smth)
 154                                  arr2s=[ft.getUndergroundDataArray()[elt[1][0]:elt[1][1]] for elt in smth1]
 155                                  arr1s.append(DataArrayDouble.Aggregate(arr2s))
 156                                  pass
 157                              pass
 158                          pass
 159                      pass
 160                  else:
 161                      for ft in fts:
 162                          smth=filter(lambda x:x[0]==NORM_ERROR,ft.getFieldSplitedByType())
 163                          arr2=DataArrayDouble.Aggregate([ft.getUndergroundDataArray()[elt[1][0][1][0]:elt[1][0][1][1]] for elt in smth])
 164                          arr1s.append(arr2)
 165                          pass
 166                      pass
 167                  arr=DataArrayDouble.Aggregate(arr1s)
 168                  if typp==ON_CELLS:
 169                      arr.renumberInPlace(o2nML[lev])
 170                  mcf=MEDCouplingFieldDouble(typp,ONE_TIME) ; mcf.setName(fieldName) ; mcf.setTime(tim,dt,it) ; mcf.setArray(arr)
 171                  mcf.setMesh(mergeMLMesh.getMeshAtLevel(lev)) ; mcf.checkCoherency()
 172                  mergeField.appendFieldNoProfileSBT(mcf)
 173                  pass
 174              pass
 175          mergeMLFields.pushField(mergeField)
 176          pass
 177      mergeMLMesh.write("merge.med",2)
 178      mergeMLFields.write("merge.med",0)
 179
 180
 181 Solution
 182 ~~~~~~~~
 183
 184 :ref:`python_testMEDLoaderSplitAndMerge1_solution`