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])