litMajeurDroite = 'inondation1_riveDroite',
litMajeurGauche = 'inondation1_riveGauche',
)
-# --- value to use for Z when the node is not in a region (used to detect problems)
+# --- 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
theStudy = salome.myStudy
theStudyId = salome.myStudyId
+import numpy as np
+
# -----------------------------------------------------------------------------
from med import medfile
fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
return <fichierMaillage>F.med : med file containing the field "BOTTOM"
"""
-
+
noms = string.split(fichierMaillage,'.')
basename = string.join(noms[:-1], '.')
fichierFMaillage = basename + 'F.med'
# --- table des numeros des noeuds
numnoe = medfile.MEDINT(nnoe)
- # --- Lecture des composantes des coordonnees des noeuds
+ # --- 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
+ # --- creation du champ
nomcha1 = "BOTTOM"
ncomp1 = 1
# --1234567890123456--
comp1 = "z "
unit1 = "m "
dtunit = ""
- medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
+ medfield.MEDfieldCr(fid, nomcha1, medfile.MED_FLOAT64,
ncomp1, comp1, unit1, dtunit, maa)
# --- ecriture du champ
medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, nnoe, valz)
# --- fermeture du fichier
medfile.MEDfileClose(fid)
+ print fichierFMaillage, " field BOTTOM OK"
# -----------------------------------------------------------------------------
-from MEDLoader import MEDLoader, MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
+import MEDLoader
+from MEDLoader import MEDCouplingFieldDouble, ON_NODES, DataArrayDouble, MEDFileMesh
def createZfield2(fichierMaillage):
"""
fichierMaillage : 2D (x,y) med file produced by SMESH and used by interpolZ.
return <fichierMaillage>F.med : med file containing the field "BOTTOM"
"""
-
+
noms = string.split(fichierMaillage,'.')
basename = string.join(noms[:-1], '.')
fichierZMaillage = basename + 'Z.med'
mm=MEDFileMesh.New(fichierZMaillage)
mm.write(fichierFMaillage,2)
MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fichierFMaillage,fieldOnNodes)
+ print fichierFMaillage, " field BOTTOM OK"
# -----------------------------------------------------------------------------
#z2 = (z - 74.0)*10
z2 = z
return z2
-
+
# -----------------------------------------------------------------------------
import SMESH
smesh = smeshBuilder.New(theStudy)
-def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef):
+def interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod = 0):
"""
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.
dicoGroupeRegion: python dictionary giving the coorespondance of mesh groups to HYDRO regions.
Key: face group name, value: region name in the HYDRO Case
zUndef: Z value to use for nodes outside the regions (there must be none if the case is correct).
+ interpolMethod: integer value, default 0 = nearest point on bathymetry, 1 = linear interpolation
return <fichierMaillage>Z.med : med file with Z value on nodes
return <fichierMaillage>F.med : an exact copy of <fichierMaillage>Z.med
"""
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 xyz
vx.append(xyz[0])
vy.append(xyz[1])
- vz = cas.GetAltitudesForPoints( vx, vy, region )
- minz = min(vz)
- maxz = max(vz)
- statz[grp.GetName()] = (minz, maxz)
-
+ 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]
