]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/atestMEDCouplingLoaderEx2.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / doc / tutorial / atestMEDCouplingLoaderEx2.rst
1
2 .. _python_testmedcouplingloaderex2_solution:
3
4 Intersection géométrique de maillages
5 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
6
7 ::
8
9         import MEDLoader as ml
10         
11         def displayVTK(m,fname):
12                 tmp = m.deepCopy()
13                 tmp.tessellate2D(0.1)
14                 tmp.writeVTK(fname)
15                 return
16
17         # Read and clean Fixe.med
18         fixe = ml.MEDFileMesh.New("Fixe.med")
19         fixm = fixe.getMeshAtLevel(0)
20         print "Nb of nodes in the file : %i " % (fixm.getNumberOfNodes())
21         fixm.mergeNodes(1e-10)
22         print "Nb of non duplicated nodes : %i" % (fixm.getNumberOfNodes())
23         # Read and clean Mobile.med
24         mobile = ml.MEDFileMesh.New("Mobile.med")
25         mobm = mobile.getMeshAtLevel(0)
26         mobm.mergeNodes(1e-10)
27         # Visualize fixm and mobm with PARAVIEW
28         fixm2 = fixm.deepCopy()        # tessellate2D() modifies the current mesh
29         fixm2.tessellate2D(0.1)
30         fixm2.writeVTK("fixm2.vtu")
31         mobm2 = mobm.deepCopy()
32         mobm2.tessellate2D(0.1)
33         mobm2.writeVTK("mobm2.vtu")
34         # mobm2 is in several pieces, take the first one
35         zonesInMobm = mobm.partitionBySpreadZone()
36         print "Nb of zones in mobm : %i" % (len(zonesInMobm))
37         zone1Mobm = mobm[zonesInMobm[0]]
38         zone1Mobm.zipCoords()
39         displayVTK(zone1Mobm, "zone1Mobm.vtu")
40         # Get cell ids from the fix part in the boudning box of zone1Mobm
41         ids2 = fixm.getCellsInBoundingBox(zone1Mobm.getBoundingBox(),1e-10)
42         partFixm = fixm[ids2]
43         partFixm.zipCoords()
44         displayVTK(partFixm,"partFixm.vtu")
45         # Intersect partFixm with zone1Mobm
46         partFixMob, iPart, iMob = ml.MEDCouplingUMesh.Intersect2DMeshes(partFixm,zone1Mobm,1e-10)
47         partFixMob.mergeNodes(1e-10)
48         # Get the part of partFixm not included in zone1Mobm using partFixMob
49         ids3 = iMob.findIdsEqual(-1)
50         partFixmWithoutZone1Mobm = partFixMob[ids3]
51         displayVTK(partFixmWithoutZone1Mobm,"partFixmWithoutZone1Mobm.vtu")
52         # Check that intersection worked properly 
53         # Check #0
54         areaPartFixm = partFixm.getMeasureField(ml.ON_CELLS).getArray()
55         areaPartFixm.abs()
56         areaPartFixMob = partFixMob.getMeasureField(ml.ON_CELLS).getArray()
57         areaPartFixMob.abs()
58         val1=areaPartFixm.accumulate()[0]
59         val2=areaPartFixMob.accumulate()[0]
60         print "Check #0 %lf == %lf with precision 1e-8? %s" % (val1,val2,str(abs(val1-val2)<1e-8))
61         # Check #1
62         areaZone1Mobm = zone1Mobm.getMeasureField(ml.ON_CELLS).getArray()
63         areaZone1Mobm.abs()
64         val3 = areaZone1Mobm.accumulate()[0]
65         ids4 = iMob.findIdsNotEqual(-1)
66         areaPartFixMob2 = areaPartFixMob[ids4]
67         val4 = areaPartFixMob2.accumulate()[0]
68         print "Check #1 %lf == %lf with precision 1e-8 ? %s" % (val3,val4,str(abs(val3-val4)<1e-8))
69         # Check #2
70         isCheck2OK = True
71         for icell in xrange(partFixm.getNumberOfCells()):
72             ids5 = iPart.findIdsEqual(icell)
73             areaOfCells = areaPartFixMob[ids5]
74             areaOfCells.abs()
75             if abs(areaOfCells.accumulate()[0] - areaPartFixm[icell]) > 1e-9:
76                 isCheck2OK = False
77                 pass
78             pass
79         print "Check #2? %s" % (str(isCheck2OK))
80         # Indicator field creation
81         f = ml.MEDCouplingFieldDouble(ml.ON_CELLS,ml.ONE_TIME)
82         m = partFixMob.deepCopy()
83         m.tessellate2D(0.1)
84         f.setMesh(m)
85         arr = ml.DataArrayDouble(partFixMob.getNumberOfCells(),1)
86         arr[iMob.findIdsEqual(-1)] = 0.
87         arr[iMob.findIdsNotEqual(-1)] = 1.
88         f.setArray(arr)
89         f.checkConsistencyLight()
90         f.setName("Zone")
91         ml.MEDCouplingFieldDouble.WriteVTK("Zone.vtu",[f])
92         # Other zones
93         zonesMobm = ml.MEDCouplingUMesh.MergeUMeshesOnSameCoords([mobm[zonesInMobm[0]], mobm[zonesInMobm[1]], mobm[zonesInMobm[5]]])
94         zonesMobm.zipCoords()
95         partFixMob2,iPart2,iMob2 = ml.MEDCouplingUMesh.Intersect2DMeshes(partFixm,zonesMobm,1e-10)
96         partFixMob2.mergeNodes(1e-10)
97         f2 = ml.MEDCouplingFieldDouble(ml.ON_CELLS, ml.ONE_TIME)
98         m2 = partFixMob2.deepCopy()
99         m2.tessellate2D(0.1)
100         f2.setMesh(m2)
101         arr = ml.DataArrayDouble(partFixMob2.getNumberOfCells(),1)
102         arr[iMob2.findIdsEqual(-1)]=0.
103         st = 0
104         end = st + len(zonesInMobm[0])
105         arr[iMob2.findIdsInRange(st,end)] = 1.
106         st += len(zonesInMobm[0]) ; 
107         end = st + len(zonesInMobm[1])
108         arr[iMob2.findIdsInRange(st,end)] = 2.
109         st += len(zonesInMobm[1])
110         end = st + len(zonesInMobm[2])
111         arr[iMob2.findIdsInRange(st,end)] = 3.
112         f2.setArray(arr)
113         f2.checkConsistencyLight()
114         f2.setName("Zone2")
115         ml.MEDCouplingFieldDouble.WriteVTK("Zone2.vtu",[f2])