]> SALOME platform Git repositories - tools/medcoupling.git/blob - doc/tutorial/medcouplingloaderex2_en.rst
Salome HOME
Merge branch 'master' of https://codev-tuleap.cea.fr/plugins/git/salome/medcoupling
[tools/medcoupling.git] / doc / tutorial / medcouplingloaderex2_en.rst
1
2 Full example 2 - RJH
3 --------------------
4
5 Two MED files are used in this case, which are (very freely) inspired by the RJH experimental reactor. 
6
7 The first file "Fixe.med" represents the 2D geometry of the static RJH without the installations.
8
9 .. image:: images/fixm.jpg
10
11 The 2nd file "Mobile.med" represent the mobile part.
12
13 .. image:: images/mobm.jpg
14
15
16 Objective
17 ~~~~~~~~~
18
19 The aim of this exercise is to use MEDCoupling to intersect those two meshes, assign a field to it and thus localize the zones.
20
21
22 Implementation start
23 ~~~~~~~~~~~~~~~~~~~~
24
25 Import the whole Python module MEDLoader (which includes MEDCoupling). ::
26
27         from MEDLoader import *
28
29 Read and repare the static mesh "Fixe.med"
30 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
31
32 With the advanced API read the whole file "Fixe.med" and call "fixm" the MEDCouplingUMEsh instance 
33 representing the static mesh. ::
34
35         fixe=MEDFileMesh.New("Fixe.med")
36         fixm=fixe.getMeshAtLevel(0)
37
38 In what follows, it is required that any two cells touching each other share the same edges.
39 As we are in nodal connectivity mode it means that common nodes have to merged. This is not the case here.
40 Merge the nodes closer than 1e-10 and assess the impact on the node count of "fixm". ::
41
42         print "nb of nodes in file : %i"%(fixm.getNumberOfNodes())
43         fixm.mergeNodes(1e-10)
44         print "nb of non duplicated nodes : %i"%(fixm.getNumberOfNodes())
45
46 Same thing for "Mobile.med" (called "mobm"). Repair it by deleting duplicated nodes. ::
47
48         mobile=MEDFileMesh.New("Mobile.med")
49         mobm=mobile.getMeshAtLevel(0)
50         mobm.mergeNodes(1e-10)
51
52
53 Repair the "mobm" mesh converting from POLYGON to QPOLYG (temporary solution)
54 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
55
56 This section will disappear in the future. 
57 The RJH mesh being more generic than TRI6 and QUAD8 we need to store cells with an intermediate type QPOLYG 
58 (Quadratic Polygon) which is the polygonal extension to the 2D cells with a dynamic edge count.
59 For now this geometrical type QPOLYG is in MEDCoupling but there is no equivalent yet in MED file (work in progress
60 at EDF).
61 The trick for now is to store QPOLYG in standard linear polygons and to convert them after reading.
62 Only "mobm" is concerned. Convert all polygonal cells in "mobm" into QPOLYG. ::
63
64         ids=mobm.giveCellsWithType(NORM_POLYGON)
65         mobm.getNodalConnectivity()[mobm.getNodalConnectivityIndex()[ids]]=NORM_QPOLYG
66         mobm.computeTypes()
67
68 Visualize "fixm" and "mobm" using ParaView. Tesselation is needed: OpenGL doesn't handle properly circle arcs 
69 and those have to be split into smaller linear segments to be able to represent them. The method MEDCouplingUMesh.tessellate2D() achieves this but modifies the mesh (non const method in C++).
70 It only take a cut fineness parameter (0.1 will suffice (angle expressed in rd)). Remember not to modify 
71 neither "fixm" nor "mobm"! ::
72
73         fixm2=fixm.deepCopy()        # tessellate2D is non const  - a mesh copy is required
74         fixm2.tessellate2D(0.1)
75         fixm2.writeVTK("fixm2.vtu")
76         mobm2=mobm.deepCopy()
77         mobm2.tessellate2D(0.1)
78         mobm2.writeVTK("mobm2.vtu")
79
80 Define a small method displayVTK() which we will use later on. ::
81
82         def displayVTK(m,fname):
83                 tmp=m.deepCopy()
84                 tmp.tessellate2D(0.1)
85                 tmp.writeVTK(fname)
86                 return
87
88 Perform reductions and identifying zones
89 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
90
91 "mobm" is made of 6 distinct parts (see image above). We only want the first part. 
92 Use MEDCouplingUMesh.partitionBySpreadZone() to partition "mobm" in zones and only 
93 extract the first zone.
94 Name this new instance "zone1Mobm", remove all orphan nodes and display. ::
95
96         zonesInMobm=mobm.partitionBySpreadZone()
97         print "number of zones in mobm : %i"%(len(zonesInMobm))
98         zone1Mobm=mobm[zonesInMobm[0]]
99         zone1Mobm.zipCoords()
100         displayVTK(zone1Mobm,"zone1Mobm.vtu")
101
102 .. image:: images/zone1Mobm.jpg
103
104 From now on we work on "zone1Mobm". We will reduce the working area of "fixm" around "zone1Mobm".
105 To achieve this: reduce "fixm" taking only "fixm" cells located in the bounding box of "zone1Mobm" (MEDCouplingUMesh.getBoundingBox() and MEDCouplingUMesh.getCellsInBoundingBox()).
106 Name this object "partFixm", remove its orphan nodes and display it. ::
107
108         ids2=fixm.getCellsInBoundingBox(zone1Mobm.getBoundingBox(),1e-10)
109         partFixm=fixm[ids2]
110         partFixm.zipCoords()
111         displayVTK(partFixm,"partFixm.vtu")
112
113 .. image:: images/partFixmAndzone1Mobm.jpg
114
115 Geometrical intersection of the two meshes
116 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
117
118 This is the core of the exercise. 
119
120 We intersect geometrically "partFixm" and "zone1Mobm". 
121 This boils down to partition in a minimal fashion "partFixm" into cells belonging either fully to 
122 "partFixm", or to "partFixm" and "zone1Mobm". Invoke the static method 
123 MEDCouplingUMesh.Intersect2DMeshes(), with "partFixm" and "zone1Mobm", and use a precision
124 of 1e-10 (merge detection threshold). 
125 This method returns 3 parameters (see API documentation) which will be called partFixMob, iPart and iMob.
126
127 In partFixMob merge common nodes with a threshold of 1e-10. ::
128
129         partFixMob,iPart,iMob=MEDCouplingUMesh.Intersect2DMeshes(partFixm,zone1Mobm,1e-10)
130         partFixMob.mergeNodes(1e-10)
131
132 Get and display partFixm part which is not in zone1Mobm. Call this mesh partFixmWithoutZone1Mobm. ::
133
134         ids3=iMob.findIdsEqual(-1)
135         partFixmWithoutZone1Mobm=partFixMob[ids3]
136         displayVTK(partFixmWithoutZone1Mobm,"partFixmWithoutZone1Mobm.vtu")
137
138 .. image:: images/partFixmWithoutZone1Mobm.jpg
139
140
141 Let's now check the result quality given by MEDCouplingUMesh.Intersect2DMeshes. 
142 Three tests will be passed:
143
144  * (check#0) the cell area sum in partFixm equals the one in partFixMob
145  * (check#1) the cell area sum in zone1Mobm equals the same sum on the cells in partFixMob whose cell ID different of -1
146  * (check#2) for each cell in partFixm, its area equals the cell area sum in partFixMob
147
148 Area is a algebraic value. The check can be performed only if all cells are correctly oriented or at least
149 all oriented consistently.
150 To check this let's inspect the areas of the 38 cells of partFixm (variable name "areaPartFixm"). ::
151
152         areaPartFixm=partFixm.getMeasureField(ON_CELLS).getArray()
153         print areaPartFixm.getValues()
154
155 All values are negative: this MED file doesn't respect the MED file convention.
156 "partFixm" being mis-oriented and the method MEDCouplingUMesh.Intersect2DMeshes() conserving the orientation, "partFixMob" is also mis-oriented.
157 To cut long story short, we perform comparison on absolute arrays. 
158 Check then that the first test check#0 is successful
159
160         areaPartFixm=partFixm.getMeasureField(ON_CELLS).getArray()
161         areaPartFixm.abs()
162         areaPartFixMob=partFixMob.getMeasureField(ON_CELLS).getArray()
163         areaPartFixMob.abs()
164         val1=areaPartFixm.accumulate()[0]
165         val2=areaPartFixMob.accumulate()[0]
166         print "Check #0 %lf == %lf a 1e-8 ? %s"%(val1,val2,str(abs(val1-val2)<1e-8))
167
168 Now check#1. Same spirit as in check#0. ::
169
170         areaZone1Mobm=zone1Mobm.getMeasureField(ON_CELLS).getArray()
171         areaZone1Mobm.abs()
172         val3=areaZone1Mobm.accumulate()[0]
173         ids4=iMob.findIdsNotEqual(-1)
174         areaPartFixMob2=areaPartFixMob[ids4]
175         val4=areaPartFixMob2.accumulate()[0]
176         print "Check #1 %lf == %lf a 1e-8 ? %s"%(val3,val4,str(abs(val3-val4)<1e-8))
177
178 Finally check#2. ::
179
180         isCheck2OK=True
181         for icell in xrange(partFixm.getNumberOfCells()):
182             ids5=iPart.findIdsEqual(icell)
183             areaOfCells=areaPartFixMob[ids5]
184             areaOfCells.abs()
185             if abs(areaOfCells.accumulate()[0]-areaPartFixm[icell])>1e-9:
186                 isCheck2OK=False
187                 pass
188             pass
189         print "Check #2? %s"%(str(isCheck2OK))
190
191 Use intersection information to create fields
192 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
193
194 We are done with partFixMob. 
195 Now create a cell field on partFixMob by setting it to 0 on the part covering only partFixm and 1 on the overlapped
196 part. Visualize it in a VTK file. ::
197
198         f=MEDCouplingFieldDouble(ON_CELLS,ONE_TIME)
199         m=partFixMob.deepCopy() ; m.tessellate2D(0.1)
200         f.setMesh(m)
201         arr=DataArrayDouble(partFixMob.getNumberOfCells(),1)
202         arr[iMob.findIdsEqual(-1)]=0.
203         arr[iMob.findIdsNotEqual(-1)]=1.
204         f.setArray(arr)
205         f.checkConsistencyLight()
206         f.setName("Zone")
207         MEDCouplingFieldDouble.WriteVTK("Zone.vtu",[f])
208
209 .. image:: images/LocationEx2.jpg
210
211 More generally take zones 0, 1 and 5. 
212 Create a cell field whose value is 0 in the zone being exclusively part of fixm,
213 1 in the zone #0, 2 in the zone #1 and 3 in the zone #5. ::
214
215         zonesMobm=MEDCouplingUMesh.MergeUMeshesOnSameCoords([mobm[zonesInMobm[0]], mobm[zonesInMobm[1]], mobm[zonesInMobm[5]]])
216         zonesMobm.zipCoords()
217         partFixMob2,iPart2,iMob2=MEDCouplingUMesh.Intersect2DMeshes(partFixm,zonesMobm,1e-10)
218         partFixMob2.mergeNodes(1e-10)
219         f2=MEDCouplingFieldDouble(ON_CELLS,ONE_TIME)
220         m2=partFixMob2.deepCopy() ; m2.tessellate2D(0.1)
221         f2.setMesh(m2)
222         arr=DataArrayDouble(partFixMob2.getNumberOfCells(),1)
223         arr[iMob2.findIdsEqual(-1)]=0.
224         st=0 ; end=st+len(zonesInMobm[0])
225         arr[iMob2.findIdsInRange(st,end)]=1.
226         st+=len(zonesInMobm[0]) ; end=st+len(zonesInMobm[1])
227         arr[iMob2.findIdsInRange(st,end)]=2.
228         st+=len(zonesInMobm[1]) ; end=st+len(zonesInMobm[2])
229         arr[iMob2.findIdsInRange(st,end)]=3.
230         f2.setArray(arr)
231         f2.checkConsistencyLight()
232         f2.setName("Zone2")
233         MEDCouplingFieldDouble.WriteVTK("Zone2.vtu",[f2])
234
235 .. image:: images/zonesMobm.jpg
236
237 Solution
238 ~~~~~~~~
239
240 :ref:`python_testmedcouplingloaderex2_solution`