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 import MEDLoader
 113 from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
 114
 115 def createZfield2(fichierMaillage):
 116   """
 117   Complete the mesh for Telemac.
 118   Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
 119   createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface.
 120   There is an alternate method based on Med file, equivalent (createZfield1).
 121   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
 122   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
 123   """
 124
 125   noms = string.split(fichierMaillage,'.')
 126   basename = string.join(noms[:-1], '.')
 127   fichierZMaillage = basename + 'Z.med'
 128   fichierFMaillage = basename + 'F.med'
 129   print fichierFMaillage
 130
 131   mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0)
 132   fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES)
 133   fieldOnNodes.setName("BOTTOM")
 134   fieldOnNodes.setMesh(mymesh)
 135   fieldOnNodes.setArray(mymesh.getCoords()[:,2])
 136   mm=MEDFileMesh.New(fichierZMaillage)
 137   mm.write(fichierFMaillage,2)
 138   MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
 139
 140 # -----------------------------------------------------------------------------
 141
 142 import HYDROPy
 143
 144 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
 145   def GetAltitudeForPoint( self, theCoordX, theCoordY ):
 146     alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
 147     z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY )    # Custom calculation takes the base value and changes it to test
 148     #z2 = (z - 74.0)*10
 149     z2 = z
 150     return z2
 151
 152 # -----------------------------------------------------------------------------
 153
 154 import SMESH
 155 import SALOMEDS
 156 from salome.smesh import smeshBuilder
 157
 158 smesh = smeshBuilder.New(theStudy)
 159
 160 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod = 0):
 161   """
 162   interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
 163   to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
 164   interpolZ must be followed by createZfield1 or createZfield2.
 165   nomCas: Calculation Case Name in module HYDRO
 166   fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
 167   dicoGroupeRegion: python dictionary giving the coorespondance of mesh groups to HYDRO regions.
 168                     Key: face group name, value: region name in the HYDRO Case
 169   zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
 170   interpolMethod: integer value, default 0 = nearest point on bathymetry, 1 = linear interpolation
 171   return <fichierMaillage>Z.med : med file with Z value on nodes
 172   return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
 173   """
 174
 175   doc = HYDROPy.HYDROData_Document.Document( theStudyId )
 176   cas = doc.FindObjectByName(nomCas)
 177   print cas
 178   custom_inter = MyInterpolator()
 179
 180   noms = string.split(fichierMaillage,'.')
 181   basename = string.join(noms[:-1], '.')
 182   fichierZMaillage = basename + 'Z.med'
 183   fichierFMaillage = basename + 'F.med'
 184   fichierFonds = basename + '.xyz'
 185
 186   print fichierMaillage
 187   print fichierZMaillage
 188   print fichierFMaillage
 189   print fichierFonds
 190
 191   regions = {}
 192   ([maillagePlat], status) = smesh.CreateMeshesFromMED(fichierMaillage)
 193   groups = maillagePlat.GetGroups()
 194
 195   grpns = [grp for grp in groups if grp.GetType() == SMESH.NODE]
 196   if len(grpns) == 0:
 197     print "Problem! There are no groups of nodes in the mesh!"
 198     print "Please create at least the groups of nodes corresponding to each region of the HYDRO case"
 199     return {}
 200
 201
 202   for grp in groups:
 203     if grp.GetType() == SMESH.NODE:
 204       grpName = grp.GetName()
 205       print grpName
 206       if grpName in dicoGroupeRegion.keys():
 207         regions[dicoGroupeRegion[grpName]] = grp
 208
 209   fo = open(fichierFonds, 'w')
 210   statz = {}
 211   for nomreg, grp in regions.iteritems():
 212     print "------- region : ", nomreg
 213     region  = doc.FindObjectByName(nomreg)
 214     #print region
 215     #region.SetInterpolator(custom_inter)
 216     nodesIds = grp.GetListOfID()
 217     #print nodesIds
 218     vx = []
 219     vy = []
 220     for nodeId in nodesIds:
 221       xyz = maillagePlat.GetNodeXYZ(nodeId)
 222       #print xyz
 223       vx.append(xyz[0])
 224       vy.append(xyz[1])
 225     vz = cas.GetAltitudesForPoints( vx, vy, region, interpolMethod )
 226     minz = min(vz)
 227     maxz = max(vz)
 228     statz[grp.GetName()] = (minz, maxz)
 229
 230
 231     for i,nodeId in enumerate(nodesIds):
 232       x = vx[i]
 233       y = vy[i]
 234       z = vz[i]
 235       if z < -9000:
 236         z = zUndef
 237       #print i, nodeId, x, y, z
 238       maillagePlat.MoveNode(nodeId, x, y, z)
 239       l = "%10.2f %10.2f %10.2f\n"%(x, y, z)
 240       fo.write(l)
 241
 242   fo.close()
 243   maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 )
 244   maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 )
 245   return statz
 246
 247 # -----------------------------------------------------------------------------
 248