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):
 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   return <fichierMaillage>Z.med : med file with Z value on nodes
 170   return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
 171   """
 172
 173   doc = HYDROPy.HYDROData_Document.Document( theStudyId )
 174   cas = doc.FindObjectByName(nomCas)
 175   print cas
 176   custom_inter = MyInterpolator()
 177
 178   noms = string.split(fichierMaillage,'.')
 179   basename = string.join(noms[:-1], '.')
 180   fichierZMaillage = basename + 'Z.med'
 181   fichierFMaillage = basename + 'F.med'
 182   fichierFonds = basename + '.xyz'
 183
 184   print fichierMaillage
 185   print fichierZMaillage
 186   print fichierFMaillage
 187   print fichierFonds
 188
 189   regions = {}
 190   ([maillagePlat], status) = smesh.CreateMeshesFromMED(fichierMaillage)
 191   groups = maillagePlat.GetGroups()
 192
 193   grpns = [grp for grp in groups if grp.GetType() == SMESH.NODE]
 194   if len(grpns) == 0:
 195     print "Problem! There are no groups of nodes in the mesh!"
 196     print "Please create at least the groups of nodes corresponding to each region of the HYDRO case"
 197     return {}
 198
 199
 200   for grp in groups:
 201     if grp.GetType() == SMESH.NODE:
 202       grpName = grp.GetName()
 203       print grpName
 204       if grpName in dicoGroupeRegion.keys():
 205         regions[dicoGroupeRegion[grpName]] = grp
 206
 207   fo = open(fichierFonds, 'w')
 208   statz = {}
 209   for nomreg, grp in regions.iteritems():
 210     print "------- region : ", nomreg
 211     region  = doc.FindObjectByName(nomreg)
 212     #print region
 213     #region.SetInterpolator(custom_inter)
 214     nodesIds = grp.GetListOfID()
 215     #print nodesIds
 216     vx = []
 217     vy = []
 218     for nodeId in nodesIds:
 219       xyz = maillagePlat.GetNodeXYZ(nodeId)
 220       #print xyz
 221       vx.append(xyz[0])
 222       vy.append(xyz[1])
 223     vz = cas.GetAltitudesForPoints( vx, vy, region )
 224     minz = min(vz)
 225     maxz = max(vz)
 226     statz[grp.GetName()] = (minz, maxz)
 227
 228
 229     for i,nodeId in enumerate(nodesIds):
 230       x = vx[i]
 231       y = vy[i]
 232       z = vz[i]
 233       if z < -9000:
 234         z = zUndef
 235       #print i, nodeId, x, y, z
 236       maillagePlat.MoveNode(nodeId, x, y, z)
 237       l = "%10.2f %10.2f %10.2f\n"%(x, y, z)
 238       fo.write(l)
 239
 240   fo.close()
 241   maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 )
 242   maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 )
 243   return statz
 244
 245 # -----------------------------------------------------------------------------
 246