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)
23 # --- add a field on nodes of type double with z values, named "BOTTOM"
24 createZfield2(fichierMaillage)
26 __revision__ = "V2.01"
27 # -----------------------------------------------------------------------------
29 # -----------------------------------------------------------------------------
34 theStudy = salome.myStudy
35 theStudyId = salome.myStudyId
38 import MEDLoader as ml
39 import MEDCoupling as mc
41 # -----------------------------------------------------------------------------
43 from med import medfile
44 from med import medmesh
45 #from med import medfilter
46 from med import medfield
47 from med import medenum
48 #from med import medprofile
49 #from med import medlocalization
50 #from med import medlink
52 def createZfield1(fichierMaillage):
54 Complete the mesh for Telemac.
55 Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
56 createZfield1 is used after interpolZ. createZfield1 is base on med file interface.
57 There is an alternate method based on MEDLoader, equivalent (createZfield2).
58 The file <fichierMaillage>F.med produced by interpolz must exist, and is modified.
59 fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
60 return <fichierMaillage>F.med : med file containing the field "BOTTOM"
63 basename = fichierMaillage[:-4]
64 fichierFMaillage = basename + 'F.med'
65 print fichierFMaillage
67 # --- ouverture du fichier
68 fid=medfile.MEDfileOpen(fichierFMaillage, medenum.MED_ACC_RDEXT)
70 maa, sdim, mdim, meshtype, desc, dtunit, sort, nstep, repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
71 print "Maillage de nom : %s , de dimension : %d , et de type %s"%(maa,mdim,meshtype)
72 print " Dimension de l'espace : %d"%(sdim)
74 # --- Combien de noeuds a lire ?
75 nnoe, chgt, trsf = medmesh.MEDmeshnEntity(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
76 medenum.MED_NODE, medenum.MED_NONE,
77 medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
80 # --- Allocations memoire
81 # --- table des coordonnees flt : (dimension * nombre de noeuds )
82 coords = medfile.MEDFLOAT(nnoe*sdim)
83 # --- table des numeros des noeuds
84 numnoe = medfile.MEDINT(nnoe)
86 # --- Lecture des composantes des coordonnees des noeuds
87 medmesh.MEDmeshNodeCoordinateRd(fid, maa, medenum.MED_NO_DT, medenum.MED_NO_IT,
88 medenum.MED_FULL_INTERLACE, coords)
89 #print "Valeur de coords : ",coords
90 valz=medfile.MEDFLOAT([z for (i,z) in enumerate(coords) if i%3==2])
91 #print "Valeur de z : ",valz
93 # --- creation du champ
96 # --1234567890123456--
100 medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
101 ncomp1, comp1, unit1, dtunit, maa)
103 # --- ecriture du champ
105 medfield.MEDfieldValueWr(fid, nomcha1, 0, medenum.MED_NO_IT, 0.0,
106 medenum.MED_NODE, medenum.MED_NONE,
107 medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
108 # --- fermeture du fichier
109 medfile.MEDfileClose(fid)
110 print fichierFMaillage, " field BOTTOM OK"
112 # -----------------------------------------------------------------------------
115 from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
117 def createZfield2(fichierMaillage):
119 Complete the mesh for Telemac.
120 Add a field on nodes, named "BOTTOM", of type double, containing z coordinates of nodes.
121 createZfield2 is used after interpolZ. createZfield1 is base on MEDLoader interface.
122 There is an alternate method based on Med file, equivalent (createZfield1).
123 fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
124 return <fichierMaillage>F.med : med file containing the field "BOTTOM"
127 basename = fichierMaillage[:-4]
128 fichierZMaillage = basename + 'Z.med'
129 fichierFMaillage = basename + 'F.med'
130 print fichierFMaillage
132 mymesh = MEDLoader.ReadUMeshFromFile(fichierZMaillage,0)
133 fieldOnNodes=MEDCouplingFieldDouble.New(ON_NODES)
134 fieldOnNodes.setName("BOTTOM")
135 fieldOnNodes.setMesh(mymesh)
136 fieldOnNodes.setArray(mymesh.getCoords()[:,2])
137 fieldOnNodes.setTime(0.0,0,-1)
138 mm=MEDFileMesh.New(fichierZMaillage)
139 mm.write(fichierFMaillage,2)
140 MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
141 print fichierFMaillage, " field BOTTOM OK"
143 # -----------------------------------------------------------------------------
147 class MyInterpolator( HYDROPy.HYDROData_IInterpolator ):
148 def GetAltitudeForPoint( self, theCoordX, theCoordY ):
149 alt_obj = HYDROPy.HYDROData_IInterpolator.GetAltitudeObject( self );
150 z = alt_obj.GetAltitudeForPoint( theCoordX, theCoordY ) # Custom calculation takes the base value and changes it to test
155 # -----------------------------------------------------------------------------
158 def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod=0, verbose=False):
160 interpolZ takes a 2D (x,y) mesh and calls the active instance of module HYDRO
161 to interpolate the bathymetry/altimetry on the mesh nodes, to produce the Z value of each node.
162 interpolZ must be followed by createZfield1 or createZfield2.
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 coorespondance of mesh groups to HYDRO regions.
166 Key: face group name, value: region name in the HYDRO Case
167 zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
168 interpolMethod: integer value, default 0 = nearest point on bathymetry, 1 = linear interpolation
169 return <fichierMaillage>Z.med : med file with Z value on nodes
170 return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
178 doc = HYDROPy.HYDROData_Document.Document(theStudyId)
179 cas = doc.FindObjectByName(nomCas)
181 custom_inter = MyInterpolator()
183 basename = fichierMaillage[:-4]
184 fichierZMaillage = basename + 'Z.med'
185 fichierFMaillage = basename + 'F.med'
186 fichierFonds = basename + '.xyz'
188 print "fichierMaillage :", fichierMaillage
189 print "fichierZMaillage :", fichierZMaillage
190 print "fichierFMaillage :", fichierFMaillage
191 print "fichierFonds :", fichierFonds
192 print "dicoGroupeRegion =", dicoGroupeRegion
194 # 1. Reads the mesh and gets the number of nodes
196 meshMEDFileRead = ml.MEDFileMesh.New(fichierMaillage)
198 # 2. Checks the names of the groups of faces
200 l_gr_faces = list(dicoGroupeRegion.keys())
201 t_groupe_n = meshMEDFileRead.getGroupsNames()
203 for gr_face_name in l_gr_faces:
204 if gr_face_name not in t_groupe_n:
205 message += "Group: '" + gr_face_name + "'\n"
208 print "Number of problems: %d" % nb_pb
212 message += "This group does"
214 message += "That %d groups do" % nb_pb
215 message += " not belongs to the mesh\n"
216 message += "Please check the names of the groups of faces corresponding to each region of the HYDRO case"
220 # 3. Get the information about the nodes
222 nbnodes = meshMEDFileRead.getNumberOfNodes()
224 print "Number of nodes: %d" % nbnodes
226 coords = meshMEDFileRead.getCoords()
231 tb_aux = np.zeros(nbnodes, dtype=np.bool)
232 bathy = np.zeros(nbnodes, dtype=np.float)
234 # 4. Exploration of every group of faces
236 for gr_face_name in l_gr_faces:
240 nomreg = dicoGroupeRegion[gr_face_name]
241 line = "------- Region: '" + nomreg + "'"
242 line += ", connected to group '" + gr_face_name + "'"
244 region = doc.FindObjectByName(nomreg)
246 #region.SetInterpolator(custom_inter)
248 # 4.2. Mesh of the group
250 mesh_of_the_group = meshMEDFileRead.getGroup(0, gr_face_name, False)
251 nbr_cells = mesh_of_the_group.getNumberOfCells()
253 print "\t. Number of cells: %d" % nbr_cells
255 # 4.3. Nodes of the meshes of the group
258 for id_elem in range(nbr_cells):
259 l_nodes = mesh_of_the_group.getNodeIdsOfCell(id_elem)
261 for id_node in l_nodes:
262 tb_aux[id_node] = True
263 np_aux = tb_aux.nonzero()
266 print "\t. Number of nodes: %d" % len(np_aux[0])
267 #print "np_aux:", np_aux
269 # 4.4. Interpolation over the nodes of the meshes of the group
273 for nodeId in np_aux[0]:
274 vx.append(coords[nodeId, 0])
275 vy.append(coords[nodeId, 1])
278 vz = cas.GetAltitudesForPoints(vx, vy, region, interpolMethod)
284 v05z = np.percentile(vz, 05)
285 v95z = np.percentile(vz, 95)
287 ligne = ".. Minimum: %f" % minz
288 ligne += ", maximum: %f" % maxz
289 ligne += ", mean: %f\n" % meanz
290 ligne += ".. stdeviation: %f" % stdz
291 ligne += ", v05z: %f" % v05z
292 ligne += ", v95z: %f" % v95z
297 statz[gr_face_name] = (minz, maxz, meanz, stdz, v05z, v95z)
298 for iaux, nodeId in enumerate(np_aux[0]):
299 bathy[nodeId] = vz[iaux]
303 #print "bathy :\n", bathy
304 np_aux_z = (bathy < -9000.).nonzero()
306 print ".. Number of nodes below the minimum: %d" % len(np_aux_z[0])
308 for iaux in np_aux_z[0]:
314 print ".. Ecriture du champ de bathymétrie sur le fichier :\n", fichierFonds
315 fo = open(fichierFonds, 'w')
316 for nodeId in range(nbnodes):
317 #maillagePlat.MoveNode(nodeId, x, y, z)
318 line = "%10.2f %10.2f %10.2f\n" % (coords[nodeId, 0], coords[nodeId, 1], bathy[nodeId])
322 # 7 Modification of the z coordinates
324 bathy_dd = mc.DataArrayDouble(np.asfarray(bathy, dtype='float'))
325 coords3D = DataArrayDouble.Meld([coords, bathy_dd])
326 #print "coords3D\n", coords3D
327 meshMEDFileRead.setCoords(coords3D)
330 print ".. Ecriture du maillage 3D sur le fichier :\n", fichierZMaillage
331 meshMEDFileRead.write(fichierZMaillage, 2)
334 print ".. Ecriture du maillage 3D avec le champ BOTTOM sur le fichier :\n", fichierFMaillage
335 meshMEDFileRead.write(fichierFMaillage, 2)