X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROTools%2FinterpolZ.py;h=318ec82979e9abe8eab84811f5a70473a42de067;hb=49d2f9dd9efdc678c33daebea5c7069ac81ba19c;hp=ec0e96c72ebca502acf71c6da5c695a5b4023e7f;hpb=e7e82da43f8a23e87044677139f83b2654da0882;p=modules%2Fhydro.git diff --git a/src/HYDROTools/interpolZ.py b/src/HYDROTools/interpolZ.py index ec0e96c7..318ec829 100755 --- a/src/HYDROTools/interpolZ.py +++ b/src/HYDROTools/interpolZ.py @@ -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' @@ -18,18 +18,13 @@ dicoGroupeRegion= dict(litMineur = 'inondation1_litMineur', 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() @@ -37,6 +32,8 @@ theStudy = salome.myStudy theStudyId = salome.myStudyId import numpy as np +import MEDLoader as ml +import MEDCoupling as mc # ----------------------------------------------------------------------------- @@ -49,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. @@ -60,8 +58,7 @@ def createZfield1(fichierMaillage): return 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 @@ -74,8 +71,8 @@ def createZfield1(fichierMaillage): # --- 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 @@ -112,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. @@ -125,8 +120,10 @@ def createZfield2(fichierMaillage): return 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 @@ -156,103 +153,244 @@ class MyInterpolator( HYDROPy.HYDROData_IInterpolator ): # ----------------------------------------------------------------------------- -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 Z.med : med file with Z value on nodes - return F.med : an exact copy of 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 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 = 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) - - statz[grp.GetName()] = (minz, maxz, meanz, stdz, v05z, v95z) - - - 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 -# ----------------------------------------------------------------------------- -