From: Paul RASCLE Date: Tue, 11 Jul 2017 07:26:29 +0000 (+0200) Subject: interpolZ sur groupes de faces, version 1 Gérald Nicolas X-Git-Tag: Salome_8_3_Hydro_1_1rc1~1^2~3 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=37a33bcfd3424c9e3ad08c3f44cbe237c1f56799;p=modules%2Fhydro.git interpolZ sur groupes de faces, version 1 Gérald Nicolas --- diff --git a/src/HYDROTools/interpolZ.py b/src/HYDROTools/interpolZ.py index ec0e96c7..ecada786 100755 --- a/src/HYDROTools/interpolZ.py +++ b/src/HYDROTools/interpolZ.py @@ -23,13 +23,11 @@ interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef) # --- add a field on nodes of type double with z values, named "BOTTOM" createZfield2(fichierMaillage) """ +__revision__ = "V2.01" # ----------------------------------------------------------------------------- -import string - # ----------------------------------------------------------------------------- -import sys import salome salome.salome_init() @@ -37,6 +35,8 @@ theStudy = salome.myStudy theStudyId = salome.myStudyId import numpy as np +import MEDLoader as ml +import MEDCoupling as mc # ----------------------------------------------------------------------------- @@ -60,8 +60,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 +73,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 @@ -125,8 +124,7 @@ def createZfield2(fichierMaillage): return F.med : med file containing the field "BOTTOM" """ - noms = string.split(fichierMaillage,'.') - basename = string.join(noms[:-1], '.') + basename = fichierMaillage[:-4] fichierZMaillage = basename + 'Z.med' fichierFMaillage = basename + 'F.med' print fichierFMaillage @@ -156,13 +154,8 @@ 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, interpolMethod=0, 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. @@ -176,83 +169,175 @@ def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod return Z.med : med file with Z value on nodes return F.med : an exact copy of Z.med """ + statz = dict() + erreur = 0 + message = "" + + while not erreur: + + doc = HYDROPy.HYDROData_Document.Document(theStudyId) + cas = doc.FindObjectByName(nomCas) + print cas + custom_inter = MyInterpolator() + + basename = fichierMaillage[:-4] + fichierZMaillage = basename + 'Z.med' + fichierFMaillage = basename + 'F.med' + fichierFonds = basename + '.xyz' + + print "fichierMaillage :", fichierMaillage + print "fichierZMaillage :", fichierZMaillage + print "fichierFMaillage :", fichierFMaillage + print "fichierFonds :", fichierFonds + print "dicoGroupeRegion =", dicoGroupeRegion +# +# 1. Reads the mesh and gets the number of nodes +# + 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 groups of faces corresponding to each region of the HYDRO case" + erreur = 2 + break +# +# 3. Get 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] +# + tb_aux = np.zeros(nbnodes, dtype=np.bool) + bathy = np.zeros(nbnodes, dtype=np.float) +# +# 4. Exploration of every group of faces +# + for gr_face_name in l_gr_faces: +# +# 4.1. Region +# + 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 +# + 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: %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 +# + statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z) + for iaux, nodeId in enumerate(np_aux[0]): + bathy[nodeId] = vz[iaux] +# +# 5. Minimum +# + #print "bathy :\n", bathy + np_aux_z = (bathy < -9000.).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 verbose: + print ".. Ecriture du champ de bathymétrie sur le fichier :\n", fichierFonds + fo = open(fichierFonds, 'w') + for nodeId in range(nbnodes): + #maillagePlat.MoveNode(nodeId, x, y, z) + line = "%10.2f %10.2f %10.2f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId]) + fo.write(line) + fo.close() +# +# 7 Modification of the z coordinates +# + bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float')) + coords3D = DataArrayDouble.Meld([coords, bathy_dd]) + #print "coords3D\n", coords3D + meshMEDFileRead.setCoords(coords3D) +# + if verbose: + print ".. Ecriture du maillage 3D sur le fichier :\n", fichierZMaillage + meshMEDFileRead.write(fichierZMaillage, 2) + + if verbose: + print ".. Ecriture du maillage 3D avec le champ BOTTOM sur le fichier :\n", fichierFMaillage + meshMEDFileRead.write(fichierFMaillage, 2) + + break + + if erreur: + print message - 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' - - 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 -# ----------------------------------------------------------------------------- -