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]                # Extract even cells
  26         m00.simplexize(0)
  27         m01=m0[1::2]
  28         m0=MEDCouplingUMesh.MergeUMeshes([m00,m01])
  29         m0.getCoords()[:]*=1/15.   # Illustrate how to quickly rescale a mesh
  30         m0.setName("mesh")
  31
  32 .. 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.
  33
  34 Create the fields "CellField" and "NodeField" at the time-stamp (5,6) corresponding to 5.6s.
  35 ::
  36
  37         CellField=MEDCouplingFieldDouble(ON_CELLS,ONE_TIME) ; CellField.setTime(5.6,5,6) ; CellField.setMesh(m0)
  38         CellField.setName("CellField")
  39         CellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))") ; CellField.getArray().setInfoOnComponent(0,"powercell [W]")
  40         NodeField=MEDCouplingFieldDouble(ON_NODES,ONE_TIME) ; NodeField.setTime(5.6,5,6) ; NodeField.setMesh(m0)
  41         NodeField.setName("NodeField")
  42         NodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))") ; NodeField.getArray().setInfoOnComponent(0,"powernode [W]")
  43
  44 "CellField" looks like this:
  45
  46 .. image:: images/SplitAndMergeCell1.jpg
  47
  48
  49 Mesh partitionning
  50 ~~~~~~~~~~~~~~~~~~
  51
  52 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()). ::
  53
  54      proc0=m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
  55      proc1=proc0.buildComplement(m0.getNumberOfCells())
  56
  57 .. image:: images/SplitAndMerge2.jpg
  58
  59 Writing into two different MED files
  60 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  61
  62 Starting with the partition above ("proc0" and "proc1") create two MED files called "proc0.med" et "proc1.med". ::
  63
  64      NodeField0=NodeField[proc0] ; CellField0=CellField[proc0] ; CellField0.setMesh(NodeField0.getMesh())
  65      NodeField1=NodeField[proc1] ; CellField1=CellField[proc1] ; CellField1.setMesh(NodeField1.getMesh())
  66
  67      proc0_fname="proc0.med"
  68      MEDLoader.WriteField(proc0_fname,NodeField0,True)
  69      MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc0_fname,CellField0)
  70
  71      proc1_fname="proc1.med"
  72      MEDLoader.WriteField(proc1_fname,NodeField1,True)
  73      MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,CellField1)
  74
  75 Reading and merging 2 MED files - Easy (but non optimal) version
  76 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  77
  78 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". ::
  79
  80      CellField0_read=MEDLoader.ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
  81      CellField1_read=MEDLoader.ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
  82      CellField_read=MEDCouplingFieldDouble.MergeFields([CellField0_read,CellField1_read])
  83
  84 .. 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 ...
  85
  86 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).
  87 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). ::
  88
  89      CellFieldCpy=CellField.deepCopy()
  90      CellFieldCpy.substractInPlaceDM(CellField_read,10,1e-12)
  91      CellFieldCpy.getArray().abs()
  92      print CellFieldCpy.getArray().isUniform(0.,1e-12)
  93
  94 Let's do the same on "NodeField". The main difference here is that redundant information is created at the boundary. ::
  95
  96      NodeField0_read=MEDLoader.ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
  97      NodeField1_read=MEDLoader.ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
  98      NodeField_read=MEDCouplingFieldDouble.MergeFields([NodeField0_read,NodeField1_read])
  99
 100 .. note:: The mesh is read a second time here, which can be damaging in terms of performance.
 101
 102 Invoke MEDCouplingUMesh.mergeNodes() on "NodeField_read" to remove duplicate nodes.
 103 Make a deep copy called  "NodeFieldCpy" from "NodeField" and call  MEDCouplingUMesh.mergeNodes(). ::
 104
 105      NodeField_read.mergeNodes(1e-10)
 106      NodeFieldCpy=NodeField.deepCopy()
 107      NodeFieldCpy.mergeNodes(1e-10)
 108
 109 .. 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.
 110
 111 Compare "NodeFieldCpy" and "NodeField_read" still using MEDCouplingFieldDouble.substractInPlaceDM(). ::
 112
 113      NodeFieldCpy.substractInPlaceDM(NodeField_read,10,1e-12)
 114      print NodeFieldCpy.getArray().isUniform(0.,1e-12)
 115
 116
 117 Read/write of two separated MED files - More complex but more efficient version
 118 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 119
 120 We show a more systematic and more general method to merge files.
 121 This is the preferred route when dealing with big files .
 122 This method adds performance but also allows to add extra information.
 123
 124 Using the advanced API read the meshes of two files "proc0.med" and "proc1.med" and aggregate the result in an MEDFileUMesh instance "mergeMLMesh".
 125 Handle all the levels (even if there is only one in the present case) using the method
 126 MEDFileUMesh.getNonEmptyLevels() on the instance coming from "proc0.med". ::
 127
 128      fileNames=["proc0.med","proc1.med"]
 129      msML=[MEDFileMesh.New(fname) for fname in fileNames]
 130      fsML=[MEDFileFields.New(fname) for fname in fileNames]
 131      mergeMLMesh=MEDFileUMesh()
 132      mergeMLFields=MEDFileFields()
 133      for lev in msML[0].getNonEmptyLevels():
 134          o2nML=len(msML[0].getNonEmptyLevels())*[None]
 135          cs=[mML.getCoords() for mML in msML]
 136          mergeMLMesh.setCoords(DataArrayDouble.Aggregate(cs))
 137          ms=[mML.getMeshAtLevel(lev) for mML in msML]
 138          m=MEDCouplingUMesh.MergeUMeshes(ms) ; m.setCoords(mergeMLMesh.getCoords())
 139          o2nML[lev]=m.sortCellsInMEDFileFrmt()
 140          mergeMLMesh.setMeshAtLevel(lev,m)
 141          pass
 142
 143      for fieldName in fsML[0].getFieldsNames():
 144          fmts=[fML[fieldName] for fML in fsML]
 145          mergeField=MEDFileFieldMultiTS()
 146          for dt,it,tim in fmts[0].getTimeSteps():
 147              fts=[fmt[dt,it] for fmt in fmts]
 148              arrs=len(fts)*[None]
 149              for typp in fts[0].getTypesOfFieldAvailable():
 150                  arr1s=[]
 151                  if typp==ON_CELLS:
 152                      for ft in fts:
 153                          for geoTyp,smth in ft.getFieldSplitedByType():
 154                              if geoTyp!=NORM_ERROR:
 155                                  smth1=filter(lambda x:x[0]==ON_CELLS,smth)
 156                                  arr2s=[ft.getUndergroundDataArray()[elt[1][0]:elt[1][1]] for elt in smth1]
 157                                  arr1s.append(DataArrayDouble.Aggregate(arr2s))
 158                                  pass
 159                              pass
 160                          pass
 161                      pass
 162                  else:
 163                      for ft in fts:
 164                          smth=filter(lambda x:x[0]==NORM_ERROR,ft.getFieldSplitedByType())
 165                          arr2=DataArrayDouble.Aggregate([ft.getUndergroundDataArray()[elt[1][0][1][0]:elt[1][0][1][1]] for elt in smth])
 166                          arr1s.append(arr2)
 167                          pass
 168                      pass
 169                  arr=DataArrayDouble.Aggregate(arr1s)
 170                  if typp==ON_CELLS:
 171                      arr.renumberInPlace(o2nML[lev])
 172                  mcf=MEDCouplingFieldDouble(typp,ONE_TIME) ; mcf.setName(fieldName) ; mcf.setTime(tim,dt,it) ; mcf.setArray(arr)
 173                  mcf.setMesh(mergeMLMesh.getMeshAtLevel(lev)) ; mcf.checkConsistencyLight()
 174                  mergeField.appendFieldNoProfileSBT(mcf)
 175                  pass
 176              pass
 177          mergeMLFields.pushField(mergeField)
 178          pass
 179      mergeMLMesh.write("merge.med",2)
 180      mergeMLFields.write("merge.med",0)
 181
 182
 183 Solution
 184 ~~~~~~~~
 185
 186 :ref:`python_testMEDLoaderSplitAndMerge1_solution`