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.deepCpy()        # tessellate2D is non const  - a mesh copy is required
  74         fixm2.tessellate2D(0.1)
  75         fixm2.writeVTK("fixm2.vtu")
  76         mobm2=mobm.deepCpy()
  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.deepCpy()
  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 achive 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.getIdsEqual(-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.getIdsNotEqual(-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.getIdsEqual(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.deepCpy() ; m.tessellate2D(0.1)
 200         f.setMesh(m)
 201         arr=DataArrayDouble(partFixMob.getNumberOfCells(),1)
 202         arr[iMob.getIdsEqual(-1)]=0.
 203         arr[iMob.getIdsNotEqual(-1)]=1.
 204         f.setArray(arr)
 205         f.checkCoherency()
 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.deepCpy() ; m2.tessellate2D(0.1)
 221         f2.setMesh(m2)
 222         arr=DataArrayDouble(partFixMob2.getNumberOfCells(),1)
 223         arr[iMob2.getIdsEqual(-1)]=0.
 224         st=0 ; end=st+len(zonesInMobm[0])
 225         arr[iMob2.getIdsInRange(st,end)]=1.
 226         st+=len(zonesInMobm[0]) ; end=st+len(zonesInMobm[1])
 227         arr[iMob2.getIdsInRange(st,end)]=2.
 228         st+=len(zonesInMobm[1]) ; end=st+len(zonesInMobm[2])
 229         arr[iMob2.getIdsInRange(st,end)]=3.
 230         f2.setArray(arr)
 231         f2.checkCoherency()
 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`