src/HYDROTools/interpolZ.py

   1 # -*- coding: utf-8 -*-
   2
   3 """
   4 example of use case of interpolZ and createZfield methods:
   5
   6 # --- case name in HYDRO
   7 nomCas = 'inondation1'
   8
   9 # --- med file 2D(x,y) of the case produced by SMESH
  10 fichierMaillage = '/home/B27118/projets/LNHE/garonne/inondation1.med'
  11
  12 # --- dictionary [med group name] = region name
  13 dicoGroupeRegion= dict(litMineur          = 'inondation1_litMineur',
  14                        litMajeurDroite    = 'inondation1_riveDroite',
  15                        litMajeurGauche    = 'inondation1_riveGauche',
  16                        )
  17 # --- value to use for Z when the node is not in a region (used to detect problems)
  18 zUndef = 90
  19
  20 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
  21 interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef)
  22
  23 # --- add a field on nodes of type double with z values, named "BOTTOM"
  24 createZfield2(fichierMaillage)
  25 """
  26 # -----------------------------------------------------------------------------
  27
  28 import string
  29
  30 # -----------------------------------------------------------------------------
  31
  32 import sys
  33 import salome
  34
  35 salome.salome_init()
  36 theStudy = salome.myStudy
  37 theStudyId = salome.myStudyId
  38
  39 # -----------------------------------------------------------------------------
  40
  41 from med import medfile
  42 from med import medmesh
  43 #from med import medfilter
  44 from med import medfield
  45 from med import medenum
  46 #from med import medprofile
  47 #from med import medlocalization
  48 #from med import medlink
  49
  50 def createZfield1(fichierMaillage):
  51   """
  52   Complete the mesh for Telemac.
  53   Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
  54   createZfield1 is used after interpolZ. createZfield1 is base on med file interface.
  55   There is an alternate method based on MEDLoader, equivalent (createZfield2).
  56   The file <fichierMaillage>F.med produced by interpolz must exist, and is modified.
  57   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
  58   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
  59   """
  60
  61   noms = string.split(fichierMaillage,'.')
  62   basename = string.join(noms[:-1], '.')
  63   fichierFMaillage = basename + 'F.med'
  64   print fichierFMaillage
  65
  66   # --- ouverture du fichier
  67   fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT)
  68
  69   maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep,  repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
  70   print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype)
  71   print "   Dimension de l'espace : %d"%(sdim)
  72
  73   # --- Combien de noeuds a lire ?
  74   nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
  75                                              medenum.MED_NODE, medenum.MED_NONE,
  76                                              medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
  77
  78   if nnoe > 0:
  79     # --- Allocations memoire
  80     # --- table des coordonnees flt : (dimension * nombre de noeuds )
  81     coords = medfile.MEDFLOAT(nnoe*sdim)
  82     # --- table des numeros des noeuds
  83     numnoe = medfile.MEDINT(nnoe)
  84
  85     # --- Lecture des composantes des coordonnees des noeuds
  86     medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
  87                                     medenum.MED_FULL_INTERLACE, coords)
  88     #print "Valeur de coords : ",coords
  89     valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
  90     #print "Valeur de z : ",valz
  91
  92     # --- creation du champ
  93     nomcha1 = "BOTTOM"
  94     ncomp1 = 1
  95           # --1234567890123456--
  96     comp1  = "z               "
  97     unit1  = "m               "
  98     dtunit = ""
  99     medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
 100                         ncomp1, comp1, unit1, dtunit, maa)
 101
 102     # --- ecriture du champ
 103
 104     medfield.MEDfieldValueWr(fid, nomcha1, medenum.MED_NO_DT, medenum.MED_NO_IT, 0.0,
 105                              medenum.MED_NODE, medenum.MED_NONE,
 106                              medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
 107   # --- fermeture du fichier
 108   medfile.MEDfileClose(fid)
 109
 110 # -----------------------------------------------------------------------------
 111
 112 from MEDLoader import MEDLoader, MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
 113
 114 def createZfield2(fichierMaillage):
 115   """
 116   Complete the mesh for Telemac.
 117   Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
 118   createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface.
 119   There is an alternate method based on Med file, equivalent (createZfield1).
 120   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
 121   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
 122   """
 123
 124   noms = string.split(fichierMaillage,'.')
 125   basename = string.join(noms[:-1], '.')
 126   fichierZMaillage = basename + 'Z.med'
 127   fichierFMaillage = basename + 'F.med'
 128   print fichierFMaillage
 129
 130   mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0)
 131   fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES)
 132   fieldOnNodes.setName("BOTTOM")
 133   fieldOnNodes.setMesh(mymesh)
 134   fieldOnNodes.setArray(mymesh.getCoords()[:,2])
 135   mm=MEDFileMesh.New(fichierZMaillage)
 136   mm.write(fichierFMaillage,2)
 137   MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
 138
 139 # -----------------------------------------------------------------------------
 140
 141 import HYDROPy
 142
 143 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
 144   def GetAltitudeForPoint( self, theCoordX, theCoordY ):
 145     alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
 146     z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY )    # Custom calculation takes the base value and changes it to test
 147     #z2 = (z - 74.0)*10
 148     z2 = z
 149     return z2
 150
 151 # -----------------------------------------------------------------------------
 152
 153 import SMESH
 154 import SALOMEDS
 155 from salome.smesh import smeshBuilder
 156
 157 smesh = smeshBuilder.New(theStudy)
 158
 159 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod = 0):
 160   """
 161   interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
 162   to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
 163   interpolZ must be followed by createZfield1 or createZfield2.
 164   nomCas: Calculation Case Name in module HYDRO
 165   fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
 166   dicoGroupeRegion: python dictionary giving the coorespondance of mesh groups to HYDRO regions.
 167                     Key: face group name, value: region name in the HYDRO Case
 168   zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
 169   interpolMethod: integer value, default 0 = nearest point on bathymetry, 1 = linear interpolation
 170   return <fichierMaillage>Z.med : med file with Z value on nodes
 171   return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
 172   """
 173
 174   doc = HYDROPy.HYDROData_Document.Document( theStudyId )
 175   cas = doc.FindObjectByName(nomCas)
 176   print cas
 177   custom_inter = MyInterpolator()
 178
 179   noms = string.split(fichierMaillage,'.')
 180   basename = string.join(noms[:-1], '.')
 181   fichierZMaillage = basename + 'Z.med'
 182   fichierFMaillage = basename + 'F.med'
 183   fichierFonds = basename + '.xyz'
 184
 185   print fichierMaillage
 186   print fichierZMaillage
 187   print fichierFMaillage
 188   print fichierFonds
 189
 190   regions = {}
 191   ([maillagePlat], status) = smesh.CreateMeshesFromMED(fichierMaillage)
 192   groups = maillagePlat.GetGroups()
 193
 194   grpns = [grp for grp in groups if grp.GetType() == SMESH.NODE]
 195   if len(grpns) == 0:
 196     print "Problem! There are no groups of nodes in the mesh!"
 197     print "Please create at least the groups of nodes corresponding to each region of the HYDRO case"
 198     return {}
 199
 200
 201   for grp in groups:
 202     if grp.GetType() == SMESH.NODE:
 203       grpName = grp.GetName()
 204       print grpName
 205       if grpName in dicoGroupeRegion.keys():
 206         regions[dicoGroupeRegion[grpName]] = grp
 207
 208   fo = open(fichierFonds, 'w')
 209   statz = {}
 210   for nomreg, grp in regions.iteritems():
 211     print "------- region : ", nomreg
 212     region  = doc.FindObjectByName(nomreg)
 213     #print region
 214     #region.SetInterpolator(custom_inter)
 215     nodesIds = grp.GetListOfID()
 216     #print nodesIds
 217     vx = []
 218     vy = []
 219     for nodeId in nodesIds:
 220       xyz = maillagePlat.GetNodeXYZ(nodeId)
 221       #print xyz
 222       vx.append(xyz[0])
 223       vy.append(xyz[1])
 224     vz = cas.GetAltitudesForPoints( vx, vy, region, interpolMethod )
 225     minz = min(vz)
 226     maxz = max(vz)
 227     statz[grp.GetName()] = (minz, maxz)
 228
 229
 230     for i,nodeId in enumerate(nodesIds):
 231       x = vx[i]
 232       y = vy[i]
 233       z = vz[i]
 234       if z < -9000:
 235         z = zUndef
 236       #print i, nodeId, x, y, z
 237       maillagePlat.MoveNode(nodeId, x, y, z)
 238       l = "%10.2f %10.2f %10.2f\n"%(x, y, z)
 239       fo.write(l)
 240
 241   fo.close()
 242   maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 )
 243   maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 )
 244   return statz
 245
 246 # -----------------------------------------------------------------------------
 247