2 Splitting and merging a MED file using MEDLoader's advanced API
3 ---------------------------------------------------------------
8 The aim of this exercise is to create a mesh with mixed geometrical type from scratch, and to associate two fields to it:
10 * a field on cells "CellField"
11 * a field on nodes "NodeField"
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)) ::
20 import medcoupling as mc
22 m0 = mc.MEDCouplingCMesh()
23 arr = mc.DataArrayDouble(31,1) ; arr.iota(0.)
25 m0 = m0.buildUnstructured()
26 m00 = m0[::2] # Extract even cells
29 m0 = mc.MEDCouplingUMesh.MergeUMeshes([m00,m01])
30 m0.getCoords()[:] *= 1/15. # Illustrate how to quickly rescale a mesh
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.
35 Create the fields "CellField" and "NodeField" at the time-stamp (5,6) corresponding to 5.6s.
39 cellField = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
40 cellField.setTime(5.6,5,6)
42 cellField.setName("CellField")
43 cellField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
44 cellField.getArray().setInfoOnComponent(0,"powercell [W]")
46 nodeField = mc.MEDCouplingFieldDouble(mc.ON_NODES,mc.ONE_TIME)
47 nodeField.setTime(5.6,5,6)
49 nodeField.setName("NodeField")
50 nodeField.fillFromAnalytic(1,"exp(-((x-1)*(x-1)+(y-1)*(y-1)))")
51 nodeField.getArray().setInfoOnComponent(0,"powernode [W]")
53 "CellField" looks like this:
55 .. image:: images/SplitAndMergeCell1.jpg
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()).
64 proc0 = m0.getCellsInBoundingBox([(0.,0.4),(0.,0.4)],1e-10)
65 proc1 = proc0.buildComplement(m0.getNumberOfCells())
67 .. image:: images/SplitAndMerge2.jpg
69 Writing into two different MED files
70 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
72 Starting with the partition above ("proc0" and "proc1") create two MED files called "proc0.med" et "proc1.med". ::
74 nodeField0 = nodeField[proc0] ; cellField0 = cellField[proc0] ; cellField0.setMesh(nodeField0.getMesh())
75 nodeField1 = nodeField[proc1] ; cellField1 = cellField[proc1] ; cellField1.setMesh(nodeField1.getMesh())
77 proc0_fname = "proc0.med"
78 mc.WriteField(proc0_fname, nodeField0, True)
79 mc.WriteFieldUsingAlreadyWrittenMesh(proc0_fname, cellField0)
81 proc1_fname = "proc1.med"
82 mc.WriteField(proc1_fname,nodeField1,True)
83 mc.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,cellField1)
85 Reading and merging 2 MED files - Easy (but non optimal) version
86 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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". ::
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])
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 ...
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). ::
99 cellFieldCpy = cellField.deepCopy()
100 cellFieldCpy.substractInPlaceDM(cellField_read,10,1e-12)
101 cellFieldCpy.getArray().abs()
102 print(cellFieldCpy.getArray().isUniform(0.,1e-12))
104 Let's do the same on "NodeField". The main difference here is that redundant information is created at the boundary. ::
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])
110 .. note:: The mesh is read a second time here, which can be damaging in terms of performance.
112 Invoke MEDCouplingUMesh.mergeNodes() on "NodeField_read" to remove duplicate nodes.
113 Make a deep copy called "NodeFieldCpy" from "NodeField" and call MEDCouplingUMesh.mergeNodes(). ::
115 nodeField_read.mergeNodes(1e-10)
116 nodeFieldCpy = nodeField.deepCopy()
117 nodeFieldCpy.mergeNodes(1e-10)
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.
121 Compare "NodeFieldCpy" and "NodeField_read" still using MEDCouplingFieldDouble.substractInPlaceDM(). ::
123 nodeFieldCpy.substractInPlaceDM(nodeField_read,10,1e-12)
124 print(nodeFieldCpy.getArray().isUniform(0.,1e-12))
127 Read/write of two separated MED files - More complex but more efficient version
128 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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.
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". ::
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)
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():
161 if typp == mc.ON_CELLS:
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))
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])
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)
187 mergeMLFields.pushField(mergeField)
189 mergeMLMesh.write("merge.med",2)
190 mergeMLFields.write("merge.med",0)
196 :ref:`python_testMEDLoaderSplitAndMerge1_solution`