1 # -*- coding: utf-8 -*-
4 example of use case of interpolZ and createZfield methods:
6 # --- case name in HYDRO
9 # --- med file 2D(x,y) of the case produced by SMESH
10 fichierMaillage = '/home/B27118/projets/LNHE/garonne/inondation1.med'
12 # --- dictionary [med group name] = region name
13 dicoGroupeRegion= dict(litMineur = 'inondation1_litMineur',
14 litMajeurDroite = 'inondation1_riveDroite',
15 litMajeurGauche = 'inondation1_riveGauche',
17 # --- value to use for Z when the node is not in a region (used to detect problems)
20 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
21 interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod)
23 __revision__ = "V2.02"
24 # -----------------------------------------------------------------------------
26 # -----------------------------------------------------------------------------
31 theStudy = salome.myStudy
32 theStudyId = salome.myStudyId
35 import MEDLoader as ml
36 import MEDCoupling as mc
38 # -----------------------------------------------------------------------------
40 from med import medfile
41 from med import medmesh
42 #from med import medfilter
43 from med import medfield
44 from med import medenum
45 #from med import medprofile
46 #from med import medlocalization
47 #from med import medlink
49 def createZfield1(fichierMaillage):
51 Complete the mesh for Telemac.
52 Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
53 createZfield1 is used after interpolZ. createZfield1 is base on med file interface.
54 There is an alternate method based on MEDLoader, equivalent (createZfield2).
55 The file <fichierMaillage>F.med produced by interpolz must exist, and is modified.
56 fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
57 return <fichierMaillage>F.med : med file containing the field "BOTTOM"
60 basename = fichierMaillage[:-4]
61 fichierFMaillage = basename + 'F.med'
62 print fichierFMaillage
64 # --- ouverture du fichier
65 fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT)
67 maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep, repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
68 print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype)
69 print " Dimension de l'espace : %d"%(sdim)
71 # --- Combien de noeuds a lire ?
72 nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
73 medenum.MED_NODE, medenum.MED_NONE,
74 medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
77 # --- Allocations memoire
78 # --- table des coordonnees flt : (dimension * nombre de noeuds )
79 coords = medfile.MEDFLOAT(nnoe*sdim)
80 # --- table des numeros des noeuds
81 numnoe = medfile.MEDINT(nnoe)
83 # --- Lecture des composantes des coordonnees des noeuds
84 medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
85 medenum.MED_FULL_INTERLACE, coords)
86 #print "Valeur de coords : ",coords
87 valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
88 #print "Valeur de z : ",valz
90 # --- creation du champ
93 # --1234567890123456--
97 medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
98 ncomp1, comp1, unit1, dtunit, maa)
100 # --- ecriture du champ
102 medfield.MEDfieldValueWr(fid, nomcha1, 0, medenum.MED_NO_IT, 0.0,
103 medenum.MED_NODE, medenum.MED_NONE,
104 medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
105 # --- fermeture du fichier
106 medfile.MEDfileClose(fid)
107 print fichierFMaillage, " field BOTTOM OK"
109 # -----------------------------------------------------------------------------
112 from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
114 def createZfield2(fichierMaillage):
116 Complete the mesh for Telemac.
117 Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
118 createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface.
119 There is an alternate method based on Med file, equivalent (createZfield1).
120 fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
121 return <fichierMaillage>F.med : med file containing the field "BOTTOM"
124 basename = fichierMaillage[:-4]
125 fichierZMaillage = basename + 'Z.med'
126 fichierFMaillage = basename + 'F.med'
127 print fichierFMaillage
129 mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0)
130 fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES)
131 fieldOnNodes.setName("BOTTOM")
132 fieldOnNodes.setMesh(mymesh)
133 fieldOnNodes.setArray(mymesh.getCoords()[:,2])
134 fieldOnNodes.setTime(0.0,0,-1)
135 mm=MEDFileMesh.New(fichierZMaillage)
136 mm.write(fichierFMaillage,2)
137 MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
138 print fichierFMaillage, " field BOTTOM OK"
140 # -----------------------------------------------------------------------------
144 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
145 def GetAltitudeForPoint( self, theCoordX, theCoordY ):
146 alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
147 z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY ) # Custom calculation takes the base value and changes it to test
152 # -----------------------------------------------------------------------------
155 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, zUndefThreshold=-9000., xyzFile=False, verbose=False):
157 interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
158 to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
161 nomCas: Calculation Case Name in module HYDRO
162 fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
163 dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions.
164 Key: face group name, value: region name in the HYDRO Case
165 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
167 interpolMethod: integer value
168 0 = nearest point on bathymetry (default)
169 1 = linear interpolation
170 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
171 zUndefThreshold: if the interpolated z is under this value, it is replaced by zUndef
172 default value is -9000.
173 xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
175 statz: statistique for z
176 Key: face group name, value: (minz, maxz, meanz, stdz, v05z, v95z)
178 return <fichierMaillage>F.med : med file with Z value on nodes and in a field "BOTTOM"
187 print "interpolMethod: %d" % interpolMethod
188 print "zUndef:", zUndef
189 print "zUndefThreshold: %f" % zUndefThreshold
191 doc = HYDROPy.HYDROData_Document.Document(theStudyId)
192 cas = doc.FindObjectByName(nomCas)
194 custom_inter = MyInterpolator()
196 basename = fichierMaillage[:-4]
197 fichierFMaillage = basename + 'F.med'
199 print "dicoGroupeRegion =", dicoGroupeRegion
200 print "fichierMaillage =", fichierMaillage
201 print "fichierFMaillage =", fichierFMaillage
203 fichierFonds = basename + '.xyz'
204 print "fichierFonds =", fichierFonds
206 # 1. Reads the mesh and gets the number of nodes
208 meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
210 # 2. Checks the names of the groups of faces
212 l_gr_faces = list(dicoGroupeRegion.keys())
213 t_groupe_n = meshMEDFileRead.getGroupsNames()
215 for gr_face_name in l_gr_faces:
216 if gr_face_name not in t_groupe_n:
217 message += "Group: '" + gr_face_name + "'\n"
220 print "Number of problems: %d" % nb_pb
224 message += "This group does"
226 message += "That %d groups do" % nb_pb
227 message += " not belongs to the mesh.\n"
228 message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case"
232 # 3. Gets the information about the nodes
234 nbnodes = meshMEDFileRead.getNumberOfNodes()
236 print "Number of nodes: %d" % nbnodes
238 coords = meshMEDFileRead.getCoords()
243 # 4. Exploration of every group of faces
245 tb_aux = np.zeros(nbnodes, dtype=np.bool)
247 bathy = np.zeros(nbnodes, dtype=np.float)
250 for gr_face_name in l_gr_faces:
252 # 4.1. Region connected to the group
254 nomreg = dicoGroupeRegion[gr_face_name]
255 line = "------- Region: '" + nomreg + "'"
256 line += ", connected to group '" + gr_face_name + "'"
258 region = doc.FindObjectByName(nomreg)
260 #region.SetInterpolator(custom_inter)
262 # 4.2. Mesh of the group
264 mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
265 nbr_cells = mesh_of_the_group.getNumberOfCells()
267 print "\t. Number of cells: %d" % nbr_cells
269 # 4.3. Nodes of the meshes of the group
270 # Every node is flagged in tb_aux
273 for id_elem in range(nbr_cells):
274 l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
276 for id_node in l_nodes:
277 tb_aux[id_node] = True
278 np_aux = tb_aux.nonzero()
281 print "\t. Number of nodes for this group: %d" % len(np_aux[0])
282 #print "np_aux:", np_aux
284 # 4.4. Interpolation over the nodes of the meshes of the group
288 for nodeId in np_aux[0]:
289 vx.append(coords[nodeId, 0])
290 vy.append(coords[nodeId, 1])
293 vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
299 v05z = np.percentile(vz, 05)
300 v95z = np.percentile(vz, 95)
302 ligne = ".. Minimum: %f" % minz
303 ligne += ", maximum: %f" % maxz
304 ligne += ", mean: %f\n" % meanz
305 ligne += ".. stdeviation: %f" % stdz
306 ligne += ", v05z: %f" % v05z
307 ligne += ", v95z: %f" % v95z
310 # 4.5. Storage of the z and of the statistics for this region
312 statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
314 for iaux, nodeId in enumerate(np_aux[0]):
315 bathy[nodeId] = vz[iaux]
318 # If the value is lower than a threshold, an "undefined" valeur is set
320 #print "bathy :\n", bathy
321 np_aux_z = (bathy < zUndefThreshold).nonzero()
323 print ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
325 for iaux in np_aux_z[0]:
333 print ".. Ecriture du champ de bathymétrie sur le fichier :\n", fichierFonds
334 fo = open(fichierFonds, 'w')
335 for nodeId in range(nbnodes):
336 line = "%10.2f %10.2f %10.2f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
341 # 7.1. Modification of the z coordinates
343 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
344 coords3D = DataArrayDouble.Meld([coords, bathy_dd])
345 #print "coords3D\n", coords3D
346 meshMEDFileRead.setCoords(coords3D)
348 # 7.2. Writes the 3D mesh
351 print ".. Ecriture du maillage 3D sur le fichier :\n", fichierFMaillage
352 meshMEDFileRead.write(fichierFMaillage, 2)
354 # 7.3. Writes the field
357 print ".. Ecriture du champ BOTTOM"
358 fieldOnNodes=ml.MEDCouplingFieldDouble(ml.ON_NODES, ml.ONE_TIME)
359 fieldOnNodes.setName("BOTTOM")
360 fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
361 fieldOnNodes.setArray(bathy_dd)
363 fMEDFile_ch_d = ml.MEDFileField1TS()
364 fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
365 fMEDFile_ch_d.write(fichierFMaillage, 0)