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