From ef44f612e7baa53c452b3ceb2b7202b8330eb234 Mon Sep 17 00:00:00 2001 From: Paul RASCLE Date: Wed, 21 Feb 2018 11:29:54 +0100 Subject: [PATCH] cleaning of interpolZ. See EDF bugtracker #16280 --- .../h011_normalCaseManualInterpolZ.py | 8 +- .../examples/h014_caseDigueManualInterpolZ.py | 8 +- .../examples/h015_normalCaseManualTelemac.py | 8 +- .../examples/h017_interpolationLineaire.py | 17 +- .../examples/h018_streamInterpolation.py | 8 +- ...h019_normalCaseManualInterpolZStrickler.py | 8 +- src/HYDROTools/interpolZ.py | 283 ++++++++---------- 7 files changed, 157 insertions(+), 183 deletions(-) mode change 100755 => 100644 src/HYDROTools/interpolZ.py diff --git a/doc/salome/examples/h011_normalCaseManualInterpolZ.py b/doc/salome/examples/h011_normalCaseManualInterpolZ.py index ffdd3b12..6f71f1a4 100644 --- a/doc/salome/examples/h011_normalCaseManualInterpolZ.py +++ b/doc/salome/examples/h011_normalCaseManualInterpolZ.py @@ -479,7 +479,7 @@ if salome.sg.hasDesktop(): # --- Z interpolation with HYDRO #---------------------- -from salome.hydrotools.interpolZ import interpolZ, createZfield2 +from salome.hydrotools.interpolZ import interpolZ from salome.hydrotools.controls import controlStatZ # --- case name in HYDRO @@ -494,9 +494,13 @@ dicoGroupeRegion= dict(litMineur = 'garonne_1_litMineur', ) # --- value to use for Z when the node is not in a region (used to detect problems) zUndef = 90 +# --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation +interpolMethod = 0 +# --- produce a 3D mesh (Z set to its value instead of 0 +m3d = True # --- Z interpolation on the bathymety/altimetry on the mesh nodes -statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef) +statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d) #print statz refstatz = {'riveDroite': (10.88, 32.61, 24.17, 5.12, 17.57, 31.33, 0.25), 'riveGauche': (7.72, 71.86, 24.51, 12.18, 12.90, 60.36, 0.25), diff --git a/doc/salome/examples/h014_caseDigueManualInterpolZ.py b/doc/salome/examples/h014_caseDigueManualInterpolZ.py index 7026406b..adfe9dcf 100644 --- a/doc/salome/examples/h014_caseDigueManualInterpolZ.py +++ b/doc/salome/examples/h014_caseDigueManualInterpolZ.py @@ -590,7 +590,7 @@ if salome.sg.hasDesktop(): # --- Z interpolation with HYDRO #---------------------- -from salome.hydrotools.interpolZ import interpolZ, createZfield2 +from salome.hydrotools.interpolZ import interpolZ from salome.hydrotools.controls import controlStatZ # --- case name in HYDRO @@ -604,9 +604,13 @@ dicoGroupeRegion= dict(litMineur = 'garonne_1_litMineur', ) # --- value to use for Z when the node is not in a region (used to detect problems) zUndef = 90 +# --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation +interpolMethod = 0 +# --- produce a 3D mesh (Z set to its value instead of 0 +m3d = True # --- Z interpolation on the bathymety/altimetry on the mesh nodes -statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef) +statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d) #print statz refstatz = {'riveDroite': (10.88, 32.61, 24.09, 5.13, 17.57, 31.33, 0.2), 'riveGauche': (7.72, 72.40, 21.59, 8.37, 16.71, 35.71, 0.2), diff --git a/doc/salome/examples/h015_normalCaseManualTelemac.py b/doc/salome/examples/h015_normalCaseManualTelemac.py index 1929b28c..93caed24 100644 --- a/doc/salome/examples/h015_normalCaseManualTelemac.py +++ b/doc/salome/examples/h015_normalCaseManualTelemac.py @@ -486,7 +486,7 @@ if salome.sg.hasDesktop(): # --- Z interpolation with HYDRO #---------------------- -from salome.hydrotools.interpolZ import interpolZ, createZfield2 +from salome.hydrotools.interpolZ import interpolZ from salome.hydrotools.controls import controlStatZ # --- case name in HYDRO @@ -502,9 +502,13 @@ dicoGroupeRegion= dict(litMineur = 'garonne_1_litMineur', ) # --- value to use for Z when the node is not in a region (used to detect problems) zUndef = 80 +# --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation +interpolMethod = 0 +# --- produce a 3D mesh (Z set to its value instead of 0 +m3d = True # --- Z interpolation on the bathymety/altimetry on the mesh nodes -statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef) +statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d) #print statz refstatz = {'riveDroite': (10.88, 32.61, 24.17, 5.12, 17.57, 31.33, 0.25), 'riveGauche': (7.72, 71.86, 24.51, 12.18, 12.90, 60.36, 0.25), diff --git a/doc/salome/examples/h017_interpolationLineaire.py b/doc/salome/examples/h017_interpolationLineaire.py index 46e4e8f6..d9be303b 100644 --- a/doc/salome/examples/h017_interpolationLineaire.py +++ b/doc/salome/examples/h017_interpolationLineaire.py @@ -342,7 +342,7 @@ if salome.sg.hasDesktop(): # --- Z interpolation with HYDRO #---------------------- -from salome.hydrotools.interpolZ import interpolZ, createZfield2 +from salome.hydrotools.interpolZ import interpolZ from salome.hydrotools.controls import controlStatZ # --- nom du cas dans HYDRO @@ -355,16 +355,15 @@ fichierMaillage = med_file dicoGroupeRegion= dict(domaine = 'etude_Reg_1', ) -# --- méthode d'interpolation sur les nuages de points de bathymétrie -# interpolMethod = 0 : interpolation au point le plus proche -# interpolMethod = 1 : interpolation linéaire de l'altitude par triangulation du nuage de points -interpolMethod = 1 - -# --- valeur de Z à prendre quand le noeud n'est pas trouvé dans la région (détection de problèmes) +# --- value to use for Z when the node is not in a region (used to detect problems) zUndef = 90 +# --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation +interpolMethod = 1 +# --- produce a 3D mesh (Z set to its value instead of 0 +m3d = True -# --- Z interpolation Z sur la bathymetrie/altimetrie aux noeuds du maillage -statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod) +# --- Z interpolation on the bathymety/altimetry on the mesh nodes +statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d) #print statz refstatz = {'domaine': (27.10, 168.28, 91.77, 46.047, 28.637, 161.17)} controlStatZ(statz, refstatz) diff --git a/doc/salome/examples/h018_streamInterpolation.py b/doc/salome/examples/h018_streamInterpolation.py index 4387d2cd..23e2d5c3 100644 --- a/doc/salome/examples/h018_streamInterpolation.py +++ b/doc/salome/examples/h018_streamInterpolation.py @@ -455,7 +455,7 @@ controlSubMeshStats(riveGauche_1, 580) # --- Z interpolation with HYDRO #---------------------- -from salome.hydrotools.interpolZ import interpolZ, createZfield2 +from salome.hydrotools.interpolZ import interpolZ from salome.hydrotools.controls import controlStatZ # --- case name in HYDRO @@ -471,9 +471,13 @@ dicoGroupeRegion= dict(litMineur = 'Reg_litMineur', ) # --- value to use for Z when the node is not in a region (used to detect problems) zUndef = 110 +# --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation +interpolMethod = 0 +# --- produce a 3D mesh (Z set to its value instead of 0 +m3d = True # --- Z interpolation on the bathymety/altimetry on the mesh nodes -statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef) +statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d) #print statz refstatz = {'riveDroite': (100.0, 100.0, 100.0, 0.0, 100.0, 100.0), 'riveGauche': (100.0, 100.0, 100.0, 0.0, 100.0, 100.0), diff --git a/doc/salome/examples/h019_normalCaseManualInterpolZStrickler.py b/doc/salome/examples/h019_normalCaseManualInterpolZStrickler.py index cad0171b..67374e70 100644 --- a/doc/salome/examples/h019_normalCaseManualInterpolZStrickler.py +++ b/doc/salome/examples/h019_normalCaseManualInterpolZStrickler.py @@ -770,7 +770,7 @@ if salome.sg.hasDesktop(): # --- Z interpolation with HYDRO #---------------------- -from salome.hydrotools.interpolZ import interpolZ, createZfield2 +from salome.hydrotools.interpolZ import interpolZ from salome.hydrotools.controls import controlStatZ # --- case name in HYDRO @@ -785,9 +785,13 @@ dicoGroupeRegion= dict(litMineur = 'garonne_1_litMineur', ) # --- value to use for Z when the node is not in a region (used to detect problems) zUndef = 90 +# --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation +interpolMethod = 0 +# --- produce a 3D mesh (Z set to its value instead of 0 +m3d = True # --- Z interpolation on the bathymety/altimetry on the mesh nodes -statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef) +statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d) #print statz refstatz = {'riveDroite': (10.88, 32.61, 24.17, 5.12, 17.57, 31.33, 0.25), 'riveGauche': (7.72, 71.86, 24.51, 12.18, 12.90, 60.36, 0.25), diff --git a/src/HYDROTools/interpolZ.py b/src/HYDROTools/interpolZ.py old mode 100755 new mode 100644 index 318ec829..9f1c978f --- a/src/HYDROTools/interpolZ.py +++ b/src/HYDROTools/interpolZ.py @@ -1,7 +1,6 @@ # -*- coding: utf-8 -*- - """ -example of use case of interpolZ: +Example of use case of interpolZ with the default values: # --- case name in HYDRO nomCas = 'inondation1' @@ -13,14 +12,11 @@ fichierMaillage = '/home/B27118/projets/LNHE/garonne/inondation1.med' dicoGroupeRegion= dict(litMineur = 'inondation1_litMineur', litMajeurDroite = 'inondation1_riveDroite', litMajeurGauche = 'inondation1_riveGauche', - ) -# --- 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, interpolMethod) +interpolZ(nomCas, fichierMaillage, dicoGroupeRegion) """ -__revision__ = "V2.04" +__revision__ = "V3.01" # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- @@ -39,113 +35,26 @@ import MEDCoupling as mc from med import medfile from med import medmesh -#from med import medfilter from med import medfield from med import medenum -#from med import medprofile -#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. - Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes. - createZfield1 is used after interpolZ. createZfield1 is base on med file interface. - There is an alternate method based on MEDLoader, equivalent (createZfield2). - The file F.med produced by interpolz must exist, and is modified. - fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ. - return F.med : med file containing the field "BOTTOM" - """ - - basename = fichierMaillage[:-4] - fichierFMaillage = basename + 'F.med' - print fichierFMaillage - - # --- ouverture du fichier - fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT) - - maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep, repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1) - print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype) - print " Dimension de l'espace : %d"%(sdim) - - # --- 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) - - if nnoe > 0: - # --- Allocations memoire - # --- table des coordonnees flt : (dimension * nombre de noeuds ) - coords = medfile.MEDFLOAT(nnoe*sdim) - # --- table des numeros des noeuds - numnoe = medfile.MEDINT(nnoe) - - # --- 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 - nomcha1 = "BOTTOM" - ncomp1 = 1 - # --1234567890123456-- - comp1 = "z " - unit1 = "m " - dtunit = "" - medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64, - ncomp1, comp1, unit1, dtunit, maa) - - # --- ecriture du champ - - 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" - -# ----------------------------------------------------------------------------- - -# La fonction createZfield2 ne sert plus à partir du 12/07/2017 -def createZfield2(fichierMaillage): - """ - Complete the mesh for Telemac. - Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes. - createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface. - There is an alternate method based on Med file, equivalent (createZfield1). - fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ. - return 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' - print fichierFMaillage - - mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0) - fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES) - 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" # ----------------------------------------------------------------------------- import HYDROPy class MyInterpolator( HYDROPy.HYDROData_IInterpolator ): + """ + Class MyInterpolator + """ + def __init__ (self) : + """ +Constructor + """ def GetAltitudeForPoint( self, theCoordX, theCoordY ): - alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self ); + """ + Function + """ + alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self ) z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY ) # Custom calculation takes the base value and changes it to test #z2 = (z - 74.0)*10 z2 = z @@ -154,7 +63,7 @@ class MyInterpolator( HYDROPy.HYDROData_IInterpolator ): # ----------------------------------------------------------------------------- -def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, xyzFile=False, verbose=False): +def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, m3d=False, 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. @@ -170,14 +79,16 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet 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). + m3d: True/False to produce a 3D mesh. Default is False. 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" + return F.med : med file with Z value in a field "BOTTOM" + Option: Z value also in Z coordinate if m3D is true + return .xyz : text file with X, Y, Z values """ statz = dict() erreur = 0 @@ -186,48 +97,61 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet while not erreur: if verbose: - print "nomCas:", nomCas - print "interpolMethod: %d" % interpolMethod - print "zUndef:", zUndef + ligne = "nomCas: %s" % nomCas + ligne += "\ninterpolMethod: %d" % interpolMethod + if (zUndef != None ): + ligne += "\nzUndef: %f" % zUndef + ligne += "\nm3d: %d" % m3d + print (ligne) doc = HYDROPy.HYDROData_Document.Document(theStudyId) cas = doc.FindObjectByName(nomCas) - print cas + print ( "cas : ", cas) custom_inter = MyInterpolator() basename = fichierMaillage[:-4] fichierFMaillage = basename + 'F.med' - print "dicoGroupeRegion =", dicoGroupeRegion - print "fichierMaillage =", fichierMaillage - print "fichierFMaillage =", fichierFMaillage + print ("dicoGroupeRegion = ", dicoGroupeRegion) + ligne = "fichierMaillage = %s" % fichierMaillage + ligne += "\nfichierFMaillage = %s" % fichierFMaillage if xyzFile: fichierFonds = basename + '.xyz' - print "fichierFonds =", fichierFonds + ligne += "\nfichierFonds = %s" % fichierFonds + print (ligne) # # 1. Reads the mesh # meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage) + if verbose: + print (meshMEDFileRead) # # 2. Checks the names of the groups of faces # - l_gr_faces = list(dicoGroupeRegion.keys()) - t_groupe_n = meshMEDFileRead.getGroupsNames() + t_group_n = meshMEDFileRead.getGroupsNames() + dicoGroupeRegion_0 = dict() nb_pb = 0 - for gr_face_name in l_gr_faces: - if gr_face_name not in t_groupe_n: + for gr_face_name in dicoGroupeRegion: + saux = gr_face_name.strip() + if saux not in t_group_n: message += "Group: '" + gr_face_name + "'\n" nb_pb += 1 + else : + dicoGroupeRegion_0[saux] = dicoGroupeRegion[gr_face_name] if verbose: - print "Number of problems: %d" % nb_pb + ligne = "Number of problems: %d" % nb_pb + print (ligne) # 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" + message += "These %d groups do" % nb_pb + message += " not belong to the mesh.\n" + message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case.\n" + message += "Groups for this file:\n" + for group_n in t_group_n : + message += "'%s'\n" % group_n erreur = 2 break # @@ -235,12 +159,26 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet # nbnodes = meshMEDFileRead.getNumberOfNodes() if verbose: - print "Number of nodes: %d" % nbnodes + ligne = "Number of nodes: %d" % nbnodes + print (ligne) # coords = meshMEDFileRead.getCoords() #print coords - #print coords[0,0] - #print coords[0,1] + if verbose: + nb_comp = coords.getNumberOfComponents() + l_info = coords.getInfoOnComponents() + ligne = "" + l_info_0=["X", "Y", "Z"] + for id_node in (0, 1, nbnodes-1): + ligne += "\nNode #%6d:" % id_node + for iaux in range(nb_comp): + if l_info[iaux]: + saux = l_info[iaux] + else: + saux = l_info_0[iaux] + ligne += " %s" % saux + ligne += "=%f" % coords[id_node, iaux] + print (ligne) # # 4. Exploration of every group of faces # @@ -249,24 +187,23 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet bathy = np.zeros(nbnodes, dtype=np.float) bathy.fill(zUndef) # - for gr_face_name in l_gr_faces: + for gr_face_name in dicoGroupeRegion_0: # # 4.1. Region connected to the group # - nomreg = dicoGroupeRegion[gr_face_name] - line = "------- Region: '" + nomreg + "'" - line += ", connected to group '" + gr_face_name + "'" - print line + nomreg = dicoGroupeRegion_0[gr_face_name] + ligne = "------- Region: '%s'" % nomreg + ligne += ", connected to group '%s'" % gr_face_name + print (ligne) 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 + ligne = "\t. Number of cells: %d" % nbr_cells + print (ligne) # # 4.3. Nodes of the meshes of the group # Every node is flagged in tb_aux @@ -280,26 +217,30 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet 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 + ligne = "\t. Number of nodes for this group: %d" % len(np_aux[0]) + print (ligne) + #print ("np_aux:", np_aux) # # 4.4. Interpolation over the nodes of the meshes of the group # - vx = [] - vy = [] + vx = list() + vy = list() for nodeId in np_aux[0]: vx.append(coords[nodeId, 0]) vy.append(coords[nodeId, 1]) - #print "vx:\n", vx - #print "vy:\n", vy + #print ("vx:\n", vx) + #print ("vy:\n", vy) +# vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod) - #print "vz:\n", vz +# + #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 @@ -307,7 +248,7 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet ligne += ".. stdeviation: %f" % stdz ligne += ", v05z: %f" % v05z ligne += ", v95z: %f" % v95z - print ligne + print (ligne) # # 4.5. Storage of the z and of the statistics for this region # @@ -320,58 +261,72 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet # 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. +# convenient "undefined" value is set. This new undefined 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# + ligne = "zUndefThreshold: %f" % zUndefThreshold + print (ligne) # #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]) + ligne = ".. Number of nodes below the minimum: %d" % len(np_aux_z[0]) + print (ligne) if len(np_aux_z[0]): for iaux in np_aux_z[0]: bathy[iaux] = zUndef # -# 6. xyz file +# 6. Option : 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() + ligne = ".. Ecriture du champ de bathymétrie sur le fichier :\n%s" % fichierFonds + print (ligne) +# + with open(fichierFonds, "w") as fo : + for nodeId in range(nbnodes): + ligne = "%11.3f %11.3f %11.3f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId]) + fo.write(ligne) # # 7. Final MED file -# 7.1. Modification of the z coordinates +# 7.1. Transformation of the bathymetry as a double array # 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 +# 7.2. If requested, modification of the z coordinate # - meshMEDFileRead.setCoords(coords3D) + if m3d: + coords3D = ml.DataArrayDouble.Meld([coords[:,0:2], 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 mesh # if verbose: - print ".. Ecriture du maillage 3D sur le fichier :\n", fichierFMaillage + if m3d: + saux = " 3D" + else: + saux = "" + ligne = ".. Ecriture du maillage" + saux + ligne += " sur le fichier :\n%s" % fichierFMaillage + print (ligne) +# meshMEDFileRead.write(fichierFMaillage, 2) # -# 7.3. Writes the field +# 7.4. Writes the field # med_field_name = "BOTTOM" if verbose: - print ".. Ecriture du champ '"+med_field_name+"'" + ligne = ".. Ecriture du champ '%s'" % med_field_name + print (ligne) +# #print "bathy_dd =\n", bathy_dd fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES) fieldOnNodes.setName(med_field_name) @@ -388,9 +343,9 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMet fMEDFile_ch_d.write(fichierFMaillage, 0) # break - +# if erreur: print message - +# return statz -- 2.39.2