# -*- coding: utf-8 -*-
+"""
+assignStrickler
+"""
+__revision__ = "V2.02"
+
import os
-import sys
-import shutil
import salome
-from med import medfile
-from med import medmesh
-from med import medfield
-from med import medenum
-
import HYDROPy
+import numpy as np
+import MEDLoader as ml
+import MEDCoupling as mc
+
# Get current study id
salome.salome_init()
theStudy = salome.myStudy
theStudyId = salome.myStudyId
-"""Bad parameters exception"""
class BadParamsError(ValueError):
+ """Bad parameters exception"""
pass
-"""
-Creates MED file with scalar field associated with the mesh of the calculation case.
-The scalar value associated with the mesh node is the Strickler coefficient of the land cover containing the node.
-case_name: calculation case name in the study
-med_file_name: path to input MED file with mesh on nodes
-output_file_name: path to output MED file with
-med_field_name: field name
-"""
-def assignStrickler(case_name, med_file_name, output_file_name, med_field_name='BOTTOM FRICTION'):
- # Check calculation case
- doc = HYDROPy.HYDROData_Document.Document( theStudyId )
- case = doc.FindObjectByName(case_name)
- if case is None:
- raise BadParamsError("Calculation case '%s' not found" % case_name)
-
- # Check input MED file
- if not os.path.exists(med_file_name):
- raise BadParamsError("Input file '%s' not exists" % med_file_name)
-
- # Copy input file to the output path, if the output path is empty - use the input path for the output
- file_path = med_file_name
- if output_file_name and (output_file_name != file_path):
- shutil.copyfile(med_file_name, output_file_name)
- file_path = output_file_name
-
- # Open input MED file
- fid = medfile.MEDfileOpen(file_path, medenum.MED_ACC_RDEXT)
-
- # Get mesh info
- mesh_name, sdim, mdim, meshtype, desc, dtunit, sort, nstep, repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
-
- nb_nodes, chgt, trsf = medmesh.MEDmeshnEntity(fid, mesh_name, medenum.MED_NO_DT, medenum.MED_NO_IT,
- medenum.MED_NODE, medenum.MED_NONE,
- medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
-
-
- # Create field
- if nb_nodes > 0:
- # Get node coordinates
- coords = medfile.MEDFLOAT(nb_nodes*sdim)
- medmesh.MEDmeshNodeCoordinateRd(fid, mesh_name, medenum.MED_NO_DT, medenum.MED_NO_IT,
- medenum.MED_FULL_INTERLACE, coords)
-
- # Get list of Strickler coefficient values for the nodes
- values = []
- x_coords = coords[0::sdim]
- y_coords = coords[1::sdim]
- #print "case.GetStricklerCoefficientForPoints"
+
+def assignStrickler(nomCas, input_file_name, output_file_name="", med_field_name='BOTTOM FRICTION', verbose=False):
+ """
+ assignStrickler creates the scalar field associated with the mesh of the calculation case that
+ represents the Strickler coefficient of the land cover.
+
+ In:
+ nomCas: Calculation Case Name in module HYDRO
+ input_file_name: med file name corresponding to the HYDRO case
+ output_file_name: med file name with the Strickler coefficient; default is the same as fichierMaillage
+ med_field_name: name of the field of the Strickler coefficient
+ default value is 'BOTTOM FRICTION'.
+ Out:
+ """
+ erreur = 0
+ message = ""
+
+ while not erreur:
+
+ if verbose:
+ print "nomCas:", nomCas
+ print "input_file_name:", input_file_name
+ print "output_file_name:", output_file_name
+ print "med_field_name:", med_field_name
+
+# 1. Controls
+# 1.1. Check calculation case
+ doc = HYDROPy.HYDROData_Document.Document(theStudyId)
+ case = doc.FindObjectByName(nomCas)
+ if case is None:
+ raise BadParamsError("Calculation case '%s' not found" % nomCas)
+
+# 1.2. Check input MED file
+ if not os.path.exists(input_file_name):
+ raise BadParamsError("Input file '%s' not exists" % input_file_name)
+
+# 2. The output file
+# 2.1. If output MED file is not defined, it is equal to the input file
+ if not output_file_name:
+ output_file_name = input_file_name
+
+# 2.2. Copy input file to the output path, if the output path is different from the input path
+ if output_file_name != input_file_name:
+ import shutil
+ shutil.copyfile(input_file_name, output_file_name)
+#
+# 3. Reads the mesh
+#
+ meshMEDFileRead = ml.MEDFileMesh.New(output_file_name)
+#
+# 4. Gets the information about the nodes
+#
+ nbnodes = meshMEDFileRead.getNumberOfNodes()
+ if verbose:
+ print "Number of nodes: %d" % nbnodes
+#
+ coords = meshMEDFileRead.getCoords()
+ #print "coords =\n", coords
+#
+# 5. Calculation of the Strickler coefficient
+ values = list()
+ x_coords = coords[:, 0]
+ y_coords = coords[:, 1]
+ #print "x_coords =\n", x_coords
+ #print "y_coords =\n", y_coords
values = case.GetStricklerCoefficientForPoints(x_coords, y_coords, 0.0, True)
-
- # Write the values to the field
- if len(values) > 0:
- values = medfile.MEDFLOAT(values)
-
- comp = "strickler coeff "
- medfield.MEDfieldCr(fid, med_field_name, medfile.MED_FLOAT64,
- 1, comp, '', '', mesh_name)
-
- medfield.MEDfieldValueWr(fid, med_field_name, 0, medenum.MED_NO_IT, 0.0,
- medenum.MED_NODE, medenum.MED_NONE,
- medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, len(values), values)
-
- print "MED field '%s' on %s nodes has been created." % (med_field_name, nb_nodes)
-
- # close MED file
- medfile.MEDfileClose(fid)
+#
+# 6. Field MED file
+#
+ coeff = mc.DataArrayDouble(np.asfarray(values, dtype='float'))
+ coeff.setInfoOnComponents(["Strickler [SI]"])
+ #print "coeff =\n", coeff
+ if verbose:
+ print ".. Ecriture du Strickler sous le nom '"+med_field_name+"'"
+ fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
+ fieldOnNodes.setName(med_field_name)
+ fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
+ fieldOnNodes.setArray(coeff)
+# 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)
+ fMEDFile_ch_d.write(output_file_name, 0)
+#
+ break
+
+ if erreur:
+ print message
+
+ return
# -*- coding: utf-8 -*-
"""
-example of use case of interpolZ and createZfield methods:
+example of use case of interpolZ:
# --- case name in HYDRO
nomCas = 'inondation1'
# --- Z interpolation on the bathymety/altimetry on the mesh nodes
interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod)
"""
-__revision__ = "V2.02"
+__revision__ = "V2.04"
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------
#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.
# -----------------------------------------------------------------------------
-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.
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'
# -----------------------------------------------------------------------------
-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.
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"
"""
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)
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)
#
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()
# 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
#
# 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)
