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