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