1 # -*- coding: utf-8 -*-
4 example of use case of interpolZ:
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.04"
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 # La fonction createZfield1 ne sert plus à partir du 12/07/2017
50 def createZfield1(fichierMaillage):
52 Complete the mesh for Telemac.
53 Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
54 createZfield1 is used after interpolZ. createZfield1 is base on med file interface.
55 There is an alternate method based on MEDLoader, equivalent (createZfield2).
56 The file <fichierMaillage>F.med produced by interpolz must exist, and is modified.
57 fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
58 return <fichierMaillage>F.med : med file containing the field "BOTTOM"
61 basename = fichierMaillage[:-4]
62 fichierFMaillage = basename + 'F.med'
63 print fichierFMaillage
65 # --- ouverture du fichier
66 fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT)
68 maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep, repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
69 print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype)
70 print " Dimension de l'espace : %d"%(sdim)
72 # --- Combien de noeuds a lire ?
73 nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
74 medenum.MED_NODE, medenum.MED_NONE,
75 medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
78 # --- Allocations memoire
79 # --- table des coordonnees flt : (dimension * nombre de noeuds )
80 coords = medfile.MEDFLOAT(nnoe*sdim)
81 # --- table des numeros des noeuds
82 numnoe = medfile.MEDINT(nnoe)
84 # --- Lecture des composantes des coordonnees des noeuds
85 medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
86 medenum.MED_FULL_INTERLACE, coords)
87 #print "Valeur de coords : ",coords
88 valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
89 #print "Valeur de z : ",valz
91 # --- creation du champ
94 # --1234567890123456--
98 medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
99 ncomp1, comp1, unit1, dtunit, maa)
101 # --- ecriture du champ
103 medfield.MEDfieldValueWr(fid, nomcha1, 0, medenum.MED_NO_IT, 0.0,
104 medenum.MED_NODE, medenum.MED_NONE,
105 medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
106 # --- fermeture du fichier
107 medfile.MEDfileClose(fid)
108 print fichierFMaillage, " field BOTTOM OK"
110 # -----------------------------------------------------------------------------
112 # La fonction createZfield2 ne sert plus à partir du 12/07/2017
113 def createZfield2(fichierMaillage):
115 Complete the mesh for Telemac.
116 Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
117 createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface.
118 There is an alternate method based on Med file, equivalent (createZfield1).
119 fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
120 return <fichierMaillage>F.med : med file containing the field "BOTTOM"
124 from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
126 basename = fichierMaillage[:-4]
127 fichierZMaillage = basename + 'Z.med'
128 fichierFMaillage = basename + 'F.med'
129 print fichierFMaillage
131 mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0)
132 fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES)
133 fieldOnNodes.setName("BOTTOM")
134 fieldOnNodes.setMesh(mymesh)
135 fieldOnNodes.setArray(mymesh.getCoords()[:,2])
136 fieldOnNodes.setTime(0.0,0,-1)
137 mm=MEDFileMesh.New(fichierZMaillage)
138 mm.write(fichierFMaillage,2)
139 MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
140 print fichierFMaillage, " field BOTTOM OK"
142 # -----------------------------------------------------------------------------
146 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
147 def GetAltitudeForPoint( self, theCoordX, theCoordY ):
148 alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
149 z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY ) # Custom calculation takes the base value and changes it to test
154 # -----------------------------------------------------------------------------
157 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef=90., interpolMethod=0, xyzFile=False, verbose=False):
159 interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
160 to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
163 nomCas: Calculation Case Name in module HYDRO
164 fichierMaillage: med file name produced by SMESH, corresponding to the HYDRO case
165 dicoGroupeRegion: python dictionary giving the correspondance of mesh groups to HYDRO regions.
167 Value: region name in the HYDRO Case
168 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
170 interpolMethod: integer value
171 0 = nearest point on bathymetry (default)
172 1 = linear interpolation
173 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
174 xyzFile: True/False to write an ascii file with xyz for every node. Default is False.
176 statz: statistique for z
178 Value: (minz, maxz, meanz, stdz, v05z, v95z)
180 return <fichierMaillage>F.med : med file with Z value on nodes and in a field "BOTTOM"
189 print "nomCas:", nomCas
190 print "interpolMethod: %d" % interpolMethod
191 print "zUndef:", zUndef
193 doc = HYDROPy.HYDROData_Document.Document(theStudyId)
194 cas = doc.FindObjectByName(nomCas)
196 custom_inter = MyInterpolator()
198 basename = fichierMaillage[:-4]
199 fichierFMaillage = basename + 'F.med'
201 print "dicoGroupeRegion =", dicoGroupeRegion
202 print "fichierMaillage =", fichierMaillage
203 print "fichierFMaillage =", fichierFMaillage
205 fichierFonds = basename + '.xyz'
206 print "fichierFonds =", fichierFonds
210 meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
212 # 2. Checks the names of the groups of faces
214 l_gr_faces = list(dicoGroupeRegion.keys())
215 t_groupe_n = meshMEDFileRead.getGroupsNames()
217 for gr_face_name in l_gr_faces:
218 if gr_face_name not in t_groupe_n:
219 message += "Group: '" + gr_face_name + "'\n"
222 print "Number of problems: %d" % nb_pb
226 message += "This group does"
228 message += "That %d groups do" % nb_pb
229 message += " not belongs to the mesh.\n"
230 message += "Please check the names of the group(s) of faces corresponding to each region of the HYDRO case"
234 # 3. Gets the information about the nodes
236 nbnodes = meshMEDFileRead.getNumberOfNodes()
238 print "Number of nodes: %d" % nbnodes
240 coords = meshMEDFileRead.getCoords()
245 # 4. Exploration of every group of faces
247 tb_aux = np.zeros(nbnodes, dtype=np.bool)
249 bathy = np.zeros(nbnodes, dtype=np.float)
252 for gr_face_name in l_gr_faces:
254 # 4.1. Region connected to the group
256 nomreg = dicoGroupeRegion[gr_face_name]
257 line = "------- Region: '" + nomreg + "'"
258 line += ", connected to group '" + gr_face_name + "'"
260 region = doc.FindObjectByName(nomreg)
262 #region.SetInterpolator(custom_inter)
264 # 4.2. Mesh of the group
266 mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
267 nbr_cells = mesh_of_the_group.getNumberOfCells()
269 print "\t. Number of cells: %d" % nbr_cells
271 # 4.3. Nodes of the meshes of the group
272 # Every node is flagged in tb_aux
275 for id_elem in range(nbr_cells):
276 l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
278 for id_node in l_nodes:
279 tb_aux[id_node] = True
280 np_aux = tb_aux.nonzero()
283 print "\t. Number of nodes for this group: %d" % len(np_aux[0])
284 #print "np_aux:", np_aux
286 # 4.4. Interpolation over the nodes of the meshes of the group
290 for nodeId in np_aux[0]:
291 vx.append(coords[nodeId, 0])
292 vy.append(coords[nodeId, 1])
295 vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
301 v05z = np.percentile(vz, 05)
302 v95z = np.percentile(vz, 95)
304 ligne = ".. Minimum: %f" % minz
305 ligne += ", maximum: %f" % maxz
306 ligne += ", mean: %f\n" % meanz
307 ligne += ".. stdeviation: %f" % stdz
308 ligne += ", v05z: %f" % v05z
309 ligne += ", v95z: %f" % v95z
312 # 4.5. Storage of the z and of the statistics for this region
314 statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
316 for iaux, nodeId in enumerate(np_aux[0]):
317 bathy[nodeId] = vz[iaux]
320 # During the interpolation, if no value is available over a node, a default value
321 # is set: -9999. It has no importance for the final computation, but if the field
322 # or the mesh is displayed, it makes huge gap. To prevent this artefact, a more
323 # convenient "undefined" value is set. This new undefinde value is given by the user.
325 # zUndefThreshold: the default value is -9000. It is tied with the value -9999. given
326 # by the interpolation when no value is defined.
328 zUndefThreshold = -9000.
330 print "zUndefThreshold: %f" % zUndefThreshold#
332 #print "bathy :\n", bathy
333 np_aux_z = (bathy < zUndefThreshold).nonzero()
335 print ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
337 for iaux in np_aux_z[0]:
345 print ".. Ecriture du champ de bathymétrie sur le fichier :\n", fichierFonds
346 fo = open(fichierFonds, 'w')
347 for nodeId in range(nbnodes):
348 line = "%10.2f %10.2f %10.2f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
353 # 7.1. Modification of the z coordinates
355 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
356 bathy_dd.setInfoOnComponents(["Z [m]"])
358 coords3D = ml.DataArrayDouble.Meld([coords, bathy_dd])
359 coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"])
360 #print "coords3D =\n", coords3D
362 meshMEDFileRead.setCoords(coords3D)
364 # 7.2. Writes the 3D mesh
367 print ".. Ecriture du maillage 3D sur le fichier :\n", fichierFMaillage
368 meshMEDFileRead.write(fichierFMaillage, 2)
370 # 7.3. Writes the field
372 med_field_name = "BOTTOM"
374 print ".. Ecriture du champ '"+med_field_name+"'"
375 #print "bathy_dd =\n", bathy_dd
376 fieldOnNodes = ml.MEDCouplingFieldDouble(ml.ON_NODES)
377 fieldOnNodes.setName(med_field_name)
378 fieldOnNodes.setMesh(meshMEDFileRead.getMeshAtLevel(0))
379 fieldOnNodes.setArray(bathy_dd)
380 # Ces valeurs d'instants sont mises pour assurer la lecture par TELEMAC
382 # numero d'itération : 0
383 # pas de numero d'ordre (-1)
384 fieldOnNodes.setTime(0.0, 0, -1)
386 fMEDFile_ch_d = ml.MEDFileField1TS()
387 fMEDFile_ch_d.setFieldNoProfileSBT(fieldOnNodes)
388 fMEDFile_ch_d.write(fichierFMaillage, 0)