Salome HOME
patch for compilation on Windows
[modules/hydro.git] / src / HYDROTools / interpolZ.py
index 94618d7217ed3a0fb93973fdd7cfa0153d1e9e4c..318ec82979e9abe8eab84811f5a70473a42de067 100755 (executable)
@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 
 """
-example of use case of interpolZ and createZfield methods:
+example of use case of interpolZ:
 
 # --- case name in HYDRO
 nomCas = 'inondation1'
@@ -20,7 +20,7 @@ zUndef = 90
 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
 interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod)
 """
-__revision__ = "V2.02"
+__revision__ = "V2.04"
 # -----------------------------------------------------------------------------
 
 # -----------------------------------------------------------------------------
@@ -46,6 +46,7 @@ from med import medenum
 #from med import medlocalization
 #from med import medlink
 
+# La fonction createZfield1 ne sert plus à partir du 12/07/2017
 def createZfield1(fichierMaillage):
   """
   Complete the mesh for Telemac.
@@ -108,9 +109,7 @@ def createZfield1(fichierMaillage):
 
 # -----------------------------------------------------------------------------
 
-import MEDLoader
-from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
-
+# La fonction createZfield2 ne sert plus à partir du 12/07/2017
 def createZfield2(fichierMaillage):
   """
   Complete the mesh for Telemac.
@@ -121,6 +120,9 @@ def createZfield2(fichierMaillage):
   return <fichierMaillage>F.med : med file containing the field "BOTTOM"
   """
 
+  import MEDLoader
+  from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
+
   basename = fichierMaillage[:-4]
   fichierZMaillage = basename + 'Z.med'
   fichierFMaillage = basename + 'F.med'
@@ -152,7 +154,7 @@ class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
 # -----------------------------------------------------------------------------
 
 
-def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, zUndefThreshold=-9000., xyzFile=False, verbose=False):
+def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, xyzFile=False, verbose=False):
   """
   interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
   to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
@@ -161,19 +163,19 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet
     nomCas: Calculation Case Name in module HYDRO
     fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
     dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions.
-                      Key: face group name, value: region name in the HYDRO Case
+                      Key: face group name
+                      Value: region name in the HYDRO Case
     zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
                      default value is 90.
     interpolMethod: integer value
                     0 = nearest point on bathymetry (default)
                     1 = linear interpolation
     zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
-    zUndefThreshold: if the interpolated z is under this value, it is replaced by zUndef
-                     default value is -9000.
     xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
   Out:
     statz: statistique for z
-                      Key: face group name, value: (minz, maxz, meanz, stdz, v05z, v95z)
+                      Key: face group name
+                      Value: (minz, maxz, meanz, stdz, v05z, v95z)
   Out:
     return <fichierMaillage>F.med : med file with Z value on nodes and in a field "BOTTOM"
   """
@@ -182,11 +184,11 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet
   message = ""
 
   while not erreur:
-    
+
     if verbose:
+      print "nomCas:", nomCas
       print "interpolMethod: %d" % interpolMethod
       print "zUndef:", zUndef
-      print "zUndefThreshold: %f" % zUndefThreshold
 
     doc = HYDROPy.HYDROData_Document.Document(theStudyId)
     cas = doc.FindObjectByName(nomCas)
@@ -203,7 +205,7 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet
       fichierFonds = basename + '.xyz'
       print "fichierFonds     =", fichierFonds
 #
-# 1. Reads the mesh and gets the number of nodes
+# 1. Reads the mesh
 #
     meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
 #
@@ -315,7 +317,17 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet
         bathy[nodeId] = vz[iaux]
 #
 # 5. Minimum:
-#    If the value is lower than a threshold, an "undefined" valeur is set
+#    During the interpolation, if no value is available over a node, a default value
+#    is set: -9999. It has no importance for the final computation, but if the field
+#    or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
+#    convenient "undefined" value is set. This new undefinde value is given by the user.
+#
+#    zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
+#                     by the interpolation when no value is defined.
+#
+    zUndefThreshold = -9000.
+    if verbose:
+      print "zUndefThreshold: %f" % zUndefThreshold#
 #
     #print "bathy :\n", bathy
     np_aux_z = (bathy < zUndefThreshold).nonzero()
@@ -341,8 +353,12 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet
 # 7.1. Modification of the z coordinates
 #
     bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
-    coords3D = DataArrayDouble.Meld([coords, bathy_dd])
-    #print "coords3D\n", coords3D
+    bathy_dd.setInfoOnComponents(["Z [m]"])
+#
+    coords3D = ml.DataArrayDouble.Meld([coords, bathy_dd])
+    coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
+    #print "coords3D =\n", coords3D
+#
     meshMEDFileRead.setCoords(coords3D)
 #
 # 7.2. Writes the 3D mesh
@@ -353,12 +369,19 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet
 #
 # 7.3. Writes the field
 #
+    med_field_name = "BOTTOM"
     if verbose:
-      print ".. Ecriture du champ BOTTOM"
-    fieldOnNodes=ml.MEDCouplingFieldDouble(ml.ON_NODES, ml.ONE_TIME)
-    fieldOnNodes.setName("BOTTOM")
+      print ".. Ecriture du champ '"+med_field_name+"'"
+    #print "bathy_dd =\n", bathy_dd
+    fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
+    fieldOnNodes.setName(med_field_name)
     fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
     fieldOnNodes.setArray(bathy_dd)
+#   Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
+#   instant = 0.0
+#   numero d'itération : 0
+#   pas de numero d'ordre (-1)
+    fieldOnNodes.setTime(0.0, 0, -1)
 #
     fMEDFile_ch_d = ml.MEDFileField1TS()
     fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)