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
 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   The file <fichierMaillage>F.med produced by interpolz must exist, and is modified.
 121   fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
 122   return <fichierMaillage>L.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   fichierLMaillage = basename + 'L.med'
 129   print fichierLMaillage
 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
 137   MEDLoader.WriteField(fichierLMaillage,fieldOnNodes,True)
 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   for grp in groups:
 194     if grp.GetType() == SMESH.NODE:
 195       grpName = grp.GetName()
 196       print grpName
 197       if grpName in dicoGroupeRegion.keys():
 198         regions[dicoGroupeRegion[grpName]] = grp
 199
 200   fo = open(fichierFonds, 'w')
 201   statz = {}
 202   for nomreg, grp in regions.iteritems():
 203     print "------- region : ", nomreg
 204     region  = doc.FindObjectByName(nomreg)
 205     #print region
 206     #region.SetInterpolator(custom_inter)
 207     nodesIds = grp.GetListOfID()
 208     #print nodesIds
 209     vx = []
 210     vy = []
 211     for nodeId in nodesIds:
 212       xyz = maillagePlat.GetNodeXYZ(nodeId)
 213       #print xyz
 214       vx.append(xyz[0])
 215       vy.append(xyz[1])
 216     vz = cas.GetAltitudesForPoints( vx, vy, region )
 217     minz = min(vz)
 218     maxz = max(vz)
 219     statz[grp.GetName()] = (minz, maxz)
 220
 221
 222     for i,nodeId in enumerate(nodesIds):
 223       x = vx[i]
 224       y = vy[i]
 225       z = vz[i]
 226       if z < -9000:
 227         z = zUndef
 228       #print i, nodeId, x, y, z
 229       maillagePlat.MoveNode(nodeId, x, y, z)
 230       l = "%10.2f %10.2f %10.2f\n"%(x, y, z)
 231       fo.write(l)
 232
 233   fo.close()
 234   maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 )
 235   maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 )
 236   return statz
 237
 238 # -----------------------------------------------------------------------------
 239