# -*- coding: utf-8 -*-
"""
-example of use case of interpolZ and createZfield methods:
+example of use case of interpolZ:
# --- case name in HYDRO
nomCas = 'inondation1'
litMajeurDroite = 'inondation1_riveDroite',
litMajeurGauche = 'inondation1_riveGauche',
)
-# --- value to use for Z when the node is not in a region (used to detect problems)
+# --- value to use for Z when the node is not in a region (used to detect problems)
zUndef = 90
# --- Z interpolation on the bathymety/altimetry on the mesh nodes
-interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef)
-
-# --- add a field on nodes of type double with z values, named "BOTTOM"
-createZfield2(fichierMaillage)
+interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod)
"""
+__revision__ = "V2.04"
# -----------------------------------------------------------------------------
-import string
-
# -----------------------------------------------------------------------------
-import sys
import salome
salome.salome_init()
theStudy = salome.myStudy
theStudyId = salome.myStudyId
+import numpy as np
+import MEDLoader as ml
+import MEDCoupling as mc
+
# -----------------------------------------------------------------------------
from med import medfile
#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.
fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
return <fichierMaillage>F.med : med file containing the field "BOTTOM"
"""
-
- noms = string.split(fichierMaillage,'.')
- basename = string.join(noms[:-1], '.')
+
+ basename = fichierMaillage[:-4]
fichierFMaillage = basename + 'F.med'
print fichierFMaillage
# --- Combien de noeuds a lire ?
nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
- medenum.MED_NODE, medenum.MED_NONE,
- medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
+ medenum.MED_NODE, medenum.MED_NONE,
+ medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
if nnoe > 0:
# --- Allocations memoire
# --- table des numeros des noeuds
numnoe = medfile.MEDINT(nnoe)
- # --- Lecture des composantes des coordonnees des noeuds
+ # --- Lecture des composantes des coordonnees des noeuds
medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
medenum.MED_FULL_INTERLACE, coords)
#print "Valeur de coords : ",coords
valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
#print "Valeur de z : ",valz
- # --- creation du champ
+ # --- creation du champ
nomcha1 = "BOTTOM"
ncomp1 = 1
# --1234567890123456--
comp1 = "z "
unit1 = "m "
dtunit = ""
- medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
+ medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
ncomp1, comp1, unit1, dtunit, maa)
# --- ecriture du champ
- medfield.MEDfieldValueWr(fid, nomcha1, medenum.MED_NO_DT, medenum.MED_NO_IT, 0.0,
+ medfield.MEDfieldValueWr(fid, nomcha1, 0, medenum.MED_NO_IT, 0.0,
medenum.MED_NODE, medenum.MED_NONE,
medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
# --- fermeture du fichier
medfile.MEDfileClose(fid)
+ print fichierFMaillage, " field BOTTOM OK"
# -----------------------------------------------------------------------------
-from MEDLoader import MEDLoader, MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
-
+# La fonction createZfield2 ne sert plus à partir du 12/07/2017
def createZfield2(fichierMaillage):
"""
Complete the mesh for Telemac.
fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
return <fichierMaillage>F.med : med file containing the field "BOTTOM"
"""
-
- noms = string.split(fichierMaillage,'.')
- basename = string.join(noms[:-1], '.')
+
+ import MEDLoader
+ from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
+
+ basename = fichierMaillage[:-4]
fichierZMaillage = basename + 'Z.med'
fichierFMaillage = basename + 'F.med'
print fichierFMaillage
fieldOnNodes.setName("BOTTOM")
fieldOnNodes.setMesh(mymesh)
fieldOnNodes.setArray(mymesh.getCoords()[:,2])
+ fieldOnNodes.setTime(0.0,0,-1)
mm=MEDFileMesh.New(fichierZMaillage)
mm.write(fichierFMaillage,2)
MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
+ print fichierFMaillage, " field BOTTOM OK"
# -----------------------------------------------------------------------------
#z2 = (z - 74.0)*10
z2 = z
return z2
-
-# -----------------------------------------------------------------------------
-import SMESH
-import SALOMEDS
-from salome.smesh import smeshBuilder
+# -----------------------------------------------------------------------------
-smesh = smeshBuilder.New(theStudy)
-def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod = 0):
+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.
- interpolZ must be followed by createZfield1 or createZfield2.
- nomCas: Calculation Case Name in module HYDRO
- fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
- dicoGroupeRegion: python dictionary giving the coorespondance of mesh groups to HYDRO regions.
- 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).
- interpolMethod: integer value, default 0 = nearest point on bathymetry, 1 = linear interpolation
- return <fichierMaillage>Z.med : med file with Z value on nodes
- return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
- """
-
- doc = HYDROPy.HYDROData_Document.Document( theStudyId )
- cas = doc.FindObjectByName(nomCas)
- print cas
- custom_inter = MyInterpolator()
- noms = string.split(fichierMaillage,'.')
- basename = string.join(noms[:-1], '.')
- fichierZMaillage = basename + 'Z.med'
- fichierFMaillage = basename + 'F.med'
- fichierFonds = basename + '.xyz'
+ In:
+ 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
+ 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).
+ 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)
+ Out:
+ return <fichierMaillage>F.med : med file with Z value on nodes and in a field "BOTTOM"
+ """
+ statz = dict()
+ erreur = 0
+ message = ""
+
+ while not erreur:
+
+ if verbose:
+ print "nomCas:", nomCas
+ print "interpolMethod: %d" % interpolMethod
+ print "zUndef:", zUndef
+
+ doc = HYDROPy.HYDROData_Document.Document(theStudyId)
+ cas = doc.FindObjectByName(nomCas)
+ print cas
+ custom_inter = MyInterpolator()
+
+ basename = fichierMaillage[:-4]
+ fichierFMaillage = basename + 'F.med'
+
+ print "dicoGroupeRegion =", dicoGroupeRegion
+ print "fichierMaillage =", fichierMaillage
+ print "fichierFMaillage =", fichierFMaillage
+ if xyzFile:
+ fichierFonds = basename + '.xyz'
+ print "fichierFonds =", fichierFonds
+#
+# 1. Reads the mesh
+#
+ meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
+#
+# 2. Checks the names of the groups of faces
+#
+ l_gr_faces = list(dicoGroupeRegion.keys())
+ t_groupe_n = meshMEDFileRead.getGroupsNames()
+ nb_pb = 0
+ for gr_face_name in l_gr_faces:
+ if gr_face_name not in t_groupe_n:
+ message += "Group: '" + gr_face_name + "'\n"
+ nb_pb += 1
+ if verbose:
+ print "Number of problems: %d" % nb_pb
+#
+ if nb_pb > 0:
+ if nb_pb == 1:
+ message += "This group does"
+ else:
+ message += "That %d groups do" % nb_pb
+ message += " not belongs to the mesh.\n"
+ message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case"
+ erreur = 2
+ break
+#
+# 3. Gets the information about the nodes
+#
+ nbnodes = meshMEDFileRead.getNumberOfNodes()
+ if verbose:
+ print "Number of nodes: %d" % nbnodes
+#
+ coords = meshMEDFileRead.getCoords()
+ #print coords
+ #print coords[0,0]
+ #print coords[0,1]
+#
+# 4. Exploration of every group of faces
+#
+ tb_aux = np.zeros(nbnodes, dtype=np.bool)
+#
+ bathy = np.zeros(nbnodes, dtype=np.float)
+ bathy.fill(zUndef)
+#
+ for gr_face_name in l_gr_faces:
+#
+# 4.1. Region connected to the group
+#
+ nomreg = dicoGroupeRegion[gr_face_name]
+ line = "------- Region: '" + nomreg + "'"
+ line += ", connected to group '" + gr_face_name + "'"
+ print line
+ region = doc.FindObjectByName(nomreg)
+ #print region
+ #region.SetInterpolator(custom_inter)
+#
+# 4.2. Mesh of the group
+#
+ mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
+ nbr_cells = mesh_of_the_group.getNumberOfCells()
+ if verbose:
+ print "\t. Number of cells: %d" % nbr_cells
+#
+# 4.3. Nodes of the meshes of the group
+# Every node is flagged in tb_aux
+#
+ tb_aux.fill(False)
+ for id_elem in range(nbr_cells):
+ l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
+ #print l_nodes
+ for id_node in l_nodes:
+ tb_aux[id_node] = True
+ np_aux = tb_aux.nonzero()
+ if len(np_aux[0]):
+ if verbose:
+ print "\t. Number of nodes for this group: %d" % len(np_aux[0])
+ #print "np_aux:", np_aux
+#
+# 4.4. Interpolation over the nodes of the meshes of the group
+#
+ vx = []
+ vy = []
+ for nodeId in np_aux[0]:
+ vx.append(coords[nodeId, 0])
+ vy.append(coords[nodeId, 1])
+ #print "vx:\n", vx
+ #print "vy:\n", vy
+ vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
+ #print "vz:\n", vz
+ minz = np.amin(vz)
+ maxz = np.amax(vz)
+ meanz = np.mean(vz)
+ stdz = np.std(vz)
+ v05z = np.percentile(vz, 05)
+ v95z = np.percentile(vz, 95)
+ if verbose:
+ ligne = ".. Minimum: %f" % minz
+ ligne += ", maximum: %f" % maxz
+ ligne += ", mean: %f\n" % meanz
+ ligne += ".. stdeviation: %f" % stdz
+ ligne += ", v05z: %f" % v05z
+ ligne += ", v95z: %f" % v95z
+ print ligne
+#
+# 4.5. Storage of the z and of the statistics for this region
+#
+ statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
+#
+ for iaux, nodeId in enumerate(np_aux[0]):
+ bathy[nodeId] = vz[iaux]
+#
+# 5. Minimum:
+# 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()
+ if verbose:
+ print ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
+ if len(np_aux_z[0]):
+ for iaux in np_aux_z[0]:
+ bathy[iaux] = zUndef
+#
+# 6. xyz file
+#
+ if xyzFile:
+#
+ if verbose:
+ print ".. Ecriture du champ de bathymétrie sur le fichier :\n", fichierFonds
+ fo = open(fichierFonds, 'w')
+ for nodeId in range(nbnodes):
+ line = "%10.2f %10.2f %10.2f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
+ fo.write(line)
+ fo.close()
+#
+# 7. Final MED file
+# 7.1. Modification of the z coordinates
+#
+ bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
+ 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
+#
+ if verbose:
+ print ".. Ecriture du maillage 3D sur le fichier :\n", fichierFMaillage
+ meshMEDFileRead.write(fichierFMaillage, 2)
+#
+# 7.3. Writes the field
+#
+ med_field_name = "BOTTOM"
+ if verbose:
+ 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)
+ fMEDFile_ch_d.write(fichierFMaillage, 0)
+#
+ break
+
+ if erreur:
+ print message
- print fichierMaillage
- print fichierZMaillage
- print fichierFMaillage
- print fichierFonds
-
- regions = {}
- ([maillagePlat], status) = smesh.CreateMeshesFromMED(fichierMaillage)
- groups = maillagePlat.GetGroups()
-
- grpns = [grp for grp in groups if grp.GetType() == SMESH.NODE]
- if len(grpns) == 0:
- print "Problem! There are no groups of nodes in the mesh!"
- print "Please create at least the groups of nodes corresponding to each region of the HYDRO case"
- return {}
-
-
- for grp in groups:
- if grp.GetType() == SMESH.NODE:
- grpName = grp.GetName()
- print grpName
- if grpName in dicoGroupeRegion.keys():
- regions[dicoGroupeRegion[grpName]] = grp
-
- fo = open(fichierFonds, 'w')
- statz = {}
- for nomreg, grp in regions.iteritems():
- print "------- region : ", nomreg
- region = doc.FindObjectByName(nomreg)
- #print region
- #region.SetInterpolator(custom_inter)
- nodesIds = grp.GetListOfID()
- #print nodesIds
- vx = []
- vy = []
- for nodeId in nodesIds:
- xyz = maillagePlat.GetNodeXYZ(nodeId)
- #print xyz
- vx.append(xyz[0])
- vy.append(xyz[1])
- vz = cas.GetAltitudesForPoints( vx, vy, region, interpolMethod )
- minz = min(vz)
- maxz = max(vz)
- statz[grp.GetName()] = (minz, maxz)
-
-
- for i,nodeId in enumerate(nodesIds):
- x = vx[i]
- y = vy[i]
- z = vz[i]
- if z < -9000:
- z = zUndef
- #print i, nodeId, x, y, z
- maillagePlat.MoveNode(nodeId, x, y, z)
- l = "%10.2f %10.2f %10.2f\n"%(x, y, z)
- fo.write(l)
-
- fo.close()
- maillagePlat.ExportMED(fichierZMaillage , 0, SMESH.MED_V2_2, 1 )
- maillagePlat.ExportMED(fichierFMaillage , 0, SMESH.MED_V2_2, 1 )
return statz
-# -----------------------------------------------------------------------------
-
