doc/salome/examples/h019_normalCaseManualInterpolZStrickler.py

   1 # -*- coding: utf-8 -*-
   2
   3 import os
   4 HYDRO_SAMPLES = os.path.join( os.environ["HYDRO_ROOT_DIR"], "bin/salome/test")
   5 HYDRO_RESOURCES = os.path.join( os.environ["HYDRO_ROOT_DIR"], "share/salome/resources/hydro")
   6 HYDRO_TEST_RESOURCES = os.path.join(os.environ["HYDRO_DIR"], "bin", "salome", "test", "tmp_test")
   7 if not os.path.isdir(HYDRO_TEST_RESOURCES):
   8     os.mkdir(HYDRO_TEST_RESOURCES)
   9
  10 import tempfile
  11 tmpdir = tempfile.mkdtemp()
  12 print("tmpdir=",tmpdir)
  13
  14 import sys
  15 import salome
  16
  17 salome.salome_init()
  18
  19 #----------------------
  20 # --- HYDRO
  21 #----------------------
  22
  23 from HYDROPy import *
  24 from PyQt5.QtCore import *
  25 from PyQt5.QtGui import *
  26
  27 from salome.hydrotools.hydroGeoMeshUtils import loadImage, GeolocaliseImageCoords, GeolocaliseImageReference
  28 from salome.hydrotools.hydroGeoMeshUtils import importPolylines, importBathymetry, createImmersibleZone, mergePolylines, getChildrenInStudy
  29 from salome.hydrotools.hydroGeoMeshUtils import loadStricklerTable, loadLandCoverMap
  30
  31 hydro_doc = HYDROData_Document.Document()
  32
  33 offsetX = 430000.
  34 offsetY = 6350000.
  35 hydro_doc.SetLocalCS( offsetX, offsetY )
  36
  37 garonne_ign_01 = loadImage(hydro_doc, os.path.join(HYDRO_SAMPLES, "garonne_ign_01.png"), "garonne_ign_01", 0)
  38 GeolocaliseImageCoords(garonne_ign_01,
  39                        QPoint(40, 817), QPoint(1325, 85),
  40                        QPointF(471562, 6.36775e+06), QPointF(489400, 6.37702e+06))
  41
  42 garonne_ign_02 = loadImage(hydro_doc, os.path.join(HYDRO_SAMPLES, "garonne_ign_02.png"), "garonne_ign_02", 1)
  43 GeolocaliseImageReference(garonne_ign_02, garonne_ign_01,
  44                           QPoint(1389, 447), QPoint(784, 481),
  45                           QPoint(631, 95), QPoint(26, 129))
  46
  47 garonnes = importPolylines(hydro_doc, os.path.join(HYDRO_SAMPLES, "garonne.shp"), True, 5)
  48 limites_domaine = importPolylines(hydro_doc, os.path.join(HYDRO_SAMPLES, "domaine.shp"), False, 6)
  49 lits_majeur = importPolylines(hydro_doc, os.path.join(HYDRO_SAMPLES, "lit_majeur.shp"), True, 7)
  50
  51 Cloud_02 = importBathymetry(hydro_doc, os.path.join(HYDRO_SAMPLES, "Cloud_02.xyz"))
  52 garonne_point_L93 = importBathymetry(hydro_doc, os.path.join(HYDRO_SAMPLES, "garonne_point_L93.xyz"))
  53
  54 domaineEtendu = createImmersibleZone(hydro_doc, "domaineEtendu", limites_domaine[0], Cloud_02, True, 2)
  55 litMajeur = createImmersibleZone(hydro_doc, "litMajeur", lits_majeur[0], garonne_point_L93, True, 3)
  56 litMineur = createImmersibleZone(hydro_doc, "litMineur", garonnes[0], garonne_point_L93, True, 4)
  57
  58 tableStrickler, codeStrickler = loadStricklerTable(hydro_doc, os.path.join(HYDRO_RESOURCES, "def_strickler_table_06.txt"))
  59 landCoverMap = loadLandCoverMap(hydro_doc, os.path.join(HYDRO_SAMPLES, 'HYDRO', 'CLC_decoupe.shp'),
  60                                 os.path.join(HYDRO_RESOURCES, "def_strickler_table_06.txt"), 10)
  61
  62 # Calculation case
  63 garonne_1 = hydro_doc.CreateObject( KIND_CALCULATION )
  64 garonne_1.SetName( "garonne_1" )
  65
  66 garonne_1.SetAssignmentMode( HYDROData_CalculationCase.MANUAL )
  67 garonne_1.AddGeometryObject( litMineur )
  68 garonne_1.AddGeometryObject( domaineEtendu )
  69 garonne_1.AddGeometryObject( litMajeur )
  70
  71 case_geom_group = domaineEtendu.GetGroup( 0 )
  72 garonne_1.AddGeometryGroup( case_geom_group )
  73 case_geom_group = litMineur.GetGroup( 0 )
  74 garonne_1.AddGeometryGroup( case_geom_group )
  75 case_geom_group = litMajeur.GetGroup( 0 )
  76 garonne_1.AddGeometryGroup( case_geom_group )
  77 garonne_1.SetBoundaryPolyline( limites_domaine[0] )
  78 garonne_1.SetStricklerTable( tableStrickler )
  79 garonne_1.SetLandCoverMap( landCoverMap )
  80
  81 # Start the algorithm of the partition and assignment
  82 garonne_1.Update()
  83 garonne_1_litMineur = hydro_doc.FindObjectByName( "garonne_1_Reg_1" )
  84 garonne_1_Zone_1 = hydro_doc.FindObjectByName( "garonne_1_Zone_1" )
  85 garonne_1_Zone_1.SetMergeType( HYDROData_Zone.Merge_ZMIN )
  86 garonne_1_Zone_1.SetColor( QColor( 192, 113, 64 ))
  87 garonne_1_litMineur.AddZone( garonne_1_Zone_1 )
  88
  89 garonne_1_riveDroite = hydro_doc.FindObjectByName( "garonne_1_Reg_2" )
  90 garonne_1_Zone_2 = hydro_doc.FindObjectByName( "garonne_1_Zone_2" )
  91 garonne_1_Zone_2.SetColor( QColor( 141, 192, 64 ))
  92 garonne_1_riveDroite.AddZone( garonne_1_Zone_2 )
  93
  94 garonne_1_Zone_3 = hydro_doc.FindObjectByName( "garonne_1_Zone_3" )
  95 garonne_1_Zone_3.SetMergeType( HYDROData_Zone.Merge_Object )
  96 Cloud_02 = hydro_doc.FindObjectByName( "Cloud_02" )
  97 garonne_1_Zone_3.SetMergeObject( Cloud_02 )
  98 garonne_1_Zone_3.SetColor( QColor( 64, 192, 77 ))
  99 garonne_1_riveDroite.AddZone( garonne_1_Zone_3 )
 100
 101 garonne_1_riveGauche = hydro_doc.FindObjectByName( "garonne_1_Reg_3" )
 102 garonne_1_Zone_4 = hydro_doc.FindObjectByName( "garonne_1_Zone_4" )
 103 garonne_1_Zone_4.SetMergeType( HYDROData_Zone.Merge_Object )
 104 Cloud_02 = hydro_doc.FindObjectByName( "Cloud_02" )
 105 garonne_1_Zone_4.SetMergeObject( Cloud_02 )
 106 garonne_1_Zone_4.SetColor( QColor( 64, 75, 192 ))
 107 garonne_1_riveGauche.AddZone( garonne_1_Zone_4 )
 108
 109 garonne_1_Zone_5 = hydro_doc.FindObjectByName( "garonne_1_Zone_5" )
 110 garonne_1_Zone_5.SetColor( QColor( 64, 192, 77 ))
 111 garonne_1_riveGauche.AddZone( garonne_1_Zone_5 )
 112
 113 garonne_1_litMineur.SetName("garonne_1_litMineur")
 114 garonne_1_riveDroite.SetName("garonne_1_riveDroite")
 115 garonne_1_riveGauche.SetName("garonne_1_riveGauche")
 116
 117 # Export of the calculation case
 118 garonne_1_entry = garonne_1.Export()
 119
 120 #----------------------
 121 # --- Geometry
 122 #----------------------
 123
 124 # Get geometry shape and print debug information
 125 import GEOM
 126 from salome.geom import geomBuilder
 127 import math
 128 import SALOMEDS
 129 from salome.hydrotools.controls import controlGeomProps
 130
 131 geompy = geomBuilder.New()
 132
 133 print("Entry:", garonne_1_entry)
 134 HYDRO_garonne_1 = salome.IDToObject( str( garonne_1_entry ) )
 135 print("Geom shape:", HYDRO_garonne_1)
 136 print("Geom shape name:", HYDRO_garonne_1.GetName())
 137
 138
 139 # --- manual definition: geometrical faces
 140
 141 children = getChildrenInStudy(HYDRO_garonne_1)
 142 garonne_riveGauche = children["garonne_1_riveGauche"]
 143 garonne_litMineur = children["garonne_1_litMineur"]
 144 garonne_riveDroite = children["garonne_1_riveDroite"]
 145
 146 controlGeomProps(geompy, garonne_riveGauche,  29149.36918,  35948828.352061)
 147 controlGeomProps(geompy, garonne_litMineur,   30337.548492,  3488480.304388)
 148 controlGeomProps(geompy, garonne_riveDroite,  32012.343241, 25998769.23615)
 149
 150 # --- manual identification of all useful edge groups (boundary conditions)
 151
 152 allEdgesIds = geompy.SubShapeAllIDs(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 153 print("allEdgesIds", allEdgesIds)
 154
 155 (isDone, ClosedFreeBoundary, OpenFreeBoundary) = geompy.GetFreeBoundary(HYDRO_garonne_1)
 156 geompy.addToStudyInFather(HYDRO_garonne_1, ClosedFreeBoundary[0], "ClosedFreeBoundary")
 157
 158 freeBoundary = geompy.ExtractShapes(ClosedFreeBoundary[0], geompy.ShapeType["EDGE"], True)
 159 freeBoundaryIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, freeBoundary[i]) for i in range(len(freeBoundary)) ]
 160 print("freeBoundaryIds", freeBoundaryIds)
 161
 162 [litMineur_droite] = geompy.GetSharedShapesMulti([garonne_riveDroite, garonne_litMineur], geompy.ShapeType["EDGE"], True)
 163 [litMineur_gauche] = geompy.GetSharedShapesMulti([garonne_riveGauche, garonne_litMineur], geompy.ShapeType["EDGE"], True)
 164 geompy.addToStudyInFather(HYDRO_garonne_1, litMineur_droite, "litMineur_droite")
 165 geompy.addToStudyInFather(HYDRO_garonne_1, litMineur_gauche, "litMineur_gauche")
 166 rives = [litMineur_droite, litMineur_gauche]
 167 rivesIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, rives[i]) for i in range(len(rives)) ]
 168 print("rivesIds", rivesIds)
 169
 170 edges_litMineur = geompy.GetSharedShapesMulti([HYDRO_garonne_1, garonne_litMineur], geompy.ShapeType["EDGE"], True)
 171 edges_riveGauche = geompy.GetSharedShapesMulti([HYDRO_garonne_1, garonne_riveGauche], geompy.ShapeType["EDGE"], True)
 172 edges_riveDroite = geompy.GetSharedShapesMulti([HYDRO_garonne_1, garonne_riveDroite], geompy.ShapeType["EDGE"], True)
 173 edges_litMineurIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_litMineur[i]) for i in range(len(edges_litMineur)) ]
 174 edges_riveGaucheIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveGauche[i]) for i in range(len(edges_riveGauche)) ]
 175 edges_riveDroiteIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveDroite[i]) for i in range(len(edges_riveDroite)) ]
 176
 177 print("edges_litMineurIds", edges_litMineurIds)
 178 print("edges_riveGaucheIds", edges_riveGaucheIds)
 179 print("edges_riveDroiteIds", edges_riveDroiteIds)
 180
 181 sectionsIds = [Id for Id in edges_litMineurIds if Id not in rivesIds]
 182 print("sectionsIds", sectionsIds)
 183 SectionsGaronne = geompy.CreateGroup(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 184 geompy.UnionIDs(SectionsGaronne, sectionsIds)
 185 geompy.addToStudyInFather(HYDRO_garonne_1, SectionsGaronne, "SectionsGaronne")
 186
 187 bordGaucheDomaineIds = [Id for Id in freeBoundaryIds if Id in edges_riveGaucheIds]
 188 bordDroiteDomaineIds = [Id for Id in freeBoundaryIds if Id in edges_riveDroiteIds]
 189 print("bordGaucheDomaineIds", bordGaucheDomaineIds)
 190 print("bordDroiteDomaineIds", bordDroiteDomaineIds)
 191 bordGaucheDomaine = geompy.CreateGroup(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 192 geompy.UnionIDs(bordGaucheDomaine, bordGaucheDomaineIds)
 193 geompy.addToStudyInFather(HYDRO_garonne_1, bordGaucheDomaine, "bordGaucheDomaine")
 194 bordDroiteDomaine = geompy.CreateGroup(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 195 geompy.UnionIDs(bordDroiteDomaine, bordDroiteDomaineIds)
 196 geompy.addToStudyInFather(HYDRO_garonne_1, bordDroiteDomaine, "bordDroiteDomaine")
 197
 198 amont = geompy.GetEdgeNearPoint(HYDRO_garonne_1, geompy.MakeVertex(46757.861314, 25833.234752, 0))
 199 aval = geompy.GetEdgeNearPoint(HYDRO_garonne_1, geompy.MakeVertex(39078.979127, 32588.627279, 0))
 200 geompy.addToStudyInFather(HYDRO_garonne_1, amont, "amont")
 201 geompy.addToStudyInFather(HYDRO_garonne_1, aval, "aval")
 202
 203 #----------------------
 204 # --- Meshing
 205 #----------------------
 206
 207 import  SMESH, SALOMEDS
 208 from salome.smesh import smeshBuilder
 209 from salome.hydrotools.controls import controlMeshStats, controlSubMeshStats
 210 import tempfile
 211
 212 smesh = smeshBuilder.New()
 213
 214 # --- algorithms and hypothesis
 215 garonne_1 = smesh.Mesh(HYDRO_garonne_1)
 216
 217 NETGEN_2D = garonne_1.Triangle(algo=smeshBuilder.NETGEN_1D2D)
 218 NETGEN_2D_Parameters = NETGEN_2D.Parameters()
 219 NETGEN_2D_Parameters.SetMaxSize( 200 )
 220 NETGEN_2D_Parameters.SetSecondOrder( 0 )
 221 NETGEN_2D_Parameters.SetOptimize( 1 )
 222 NETGEN_2D_Parameters.SetFineness( 4 )
 223 NETGEN_2D_Parameters.SetMinSize( 50 )
 224 NETGEN_2D_Parameters.SetUseSurfaceCurvature( 1 )
 225 NETGEN_2D_Parameters.SetFuseEdges( 1 )
 226 NETGEN_2D_Parameters.SetQuadAllowed( 0 )
 227
 228 algo2D_litMineur = garonne_1.Quadrangle(algo=smeshBuilder.QUAD_MA_PROJ,geom=garonne_litMineur)
 229 algo1D_litMineur = garonne_1.Segment(geom=garonne_litMineur)
 230 hypo1D_litMineur = algo1D_litMineur.LocalLength(100,None,1e-07)
 231 subMesh_litMineur = algo1D_litMineur.GetSubMesh()
 232 smesh.SetName(subMesh_litMineur, "litMineur")
 233
 234 algo1D_SectionsGaronne = garonne_1.Segment(geom=SectionsGaronne)
 235 hypo1D_SectionsGaronne = algo1D_SectionsGaronne.NumberOfSegments(8)
 236 hypo1D_SectionsGaronne.SetDistrType( 0 )
 237 subMesh_SectionsGaronne = algo1D_SectionsGaronne.GetSubMesh()
 238 smesh.SetName(subMesh_SectionsGaronne, "SectionsGaronne")
 239
 240 isDone = garonne_1.SetMeshOrder( [ [ subMesh_SectionsGaronne, subMesh_litMineur ] ])
 241
 242 # --- compute mesh
 243 isDone = garonne_1.Compute()
 244 isDone = garonne_1.SplitQuadObject( garonne_1, 1 )
 245 isDone = garonne_1.ReorientObject( garonne_1 )
 246
 247 # --- geometrical groups of faces
 248 riveGauche_1 = garonne_1.GroupOnGeom(garonne_riveGauche,'riveGauche',SMESH.FACE)
 249 litMineur_1 = garonne_1.GroupOnGeom(garonne_litMineur,'litMineur',SMESH.FACE)
 250 riveDroite_1 = garonne_1.GroupOnGeom(garonne_riveDroite,'riveDroite',SMESH.FACE)
 251
 252 # --- geometrical groups of edges
 253
 254 ClosedFreeBoundary_1 = garonne_1.GroupOnGeom(ClosedFreeBoundary[0],'ClosedFreeBoundary',SMESH.EDGE)
 255 litMineur_droite_1 = garonne_1.GroupOnGeom(litMineur_droite,'litMineur_droite',SMESH.EDGE)
 256 litMineur_gauche_1 = garonne_1.GroupOnGeom(litMineur_gauche,'litMineur_gauche',SMESH.EDGE)
 257 SectionsGaronne_1 = garonne_1.GroupOnGeom(SectionsGaronne,'SectionsGaronne',SMESH.EDGE)
 258 bordGaucheDomaine_1 = garonne_1.GroupOnGeom(bordGaucheDomaine,'bordGaucheDomaine',SMESH.EDGE)
 259 bordDroiteDomaine_1 = garonne_1.GroupOnGeom(bordDroiteDomaine,'bordDroiteDomaine',SMESH.EDGE)
 260 amont_1 = garonne_1.GroupOnGeom(amont,'amont',SMESH.EDGE)
 261 aval_1 = garonne_1.GroupOnGeom(aval,'aval',SMESH.EDGE)
 262
 263 # --- geometrical groups of nodes
 264
 265 garonne_1_litMineur_2 = garonne_1.GroupOnGeom(garonne_litMineur,'garonne_1_litMineur',SMESH.NODE)
 266 garonne_1_riveDroite_2 = garonne_1.GroupOnGeom(garonne_riveDroite,'garonne_1_riveDroite',SMESH.NODE)
 267 garonne_1_riveGauche_2 = garonne_1.GroupOnGeom(garonne_riveGauche,'garonne_1_riveGauche',SMESH.NODE)
 268 ClosedFreeBoundary_2 = garonne_1.GroupOnGeom(ClosedFreeBoundary[0],'ClosedFreeBoundary',SMESH.NODE)
 269 litMineur_droite_2 = garonne_1.GroupOnGeom(litMineur_droite,'litMineur_droite',SMESH.NODE)
 270 litMineur_gauche_2 = garonne_1.GroupOnGeom(litMineur_gauche,'litMineur_gauche',SMESH.NODE)
 271 SectionsGaronne_2 = garonne_1.GroupOnGeom(SectionsGaronne,'SectionsGaronne',SMESH.NODE)
 272 bordGaucheDomaine_2 = garonne_1.GroupOnGeom(bordGaucheDomaine,'bordGaucheDomaine',SMESH.NODE)
 273 bordDroiteDomaine_2 = garonne_1.GroupOnGeom(bordDroiteDomaine,'bordDroiteDomaine',SMESH.NODE)
 274 amont_2 = garonne_1.GroupOnGeom(amont,'amont',SMESH.NODE)
 275 aval_2 = garonne_1.GroupOnGeom(aval,'aval',SMESH.NODE)
 276
 277 garonne_1.SetAutoColor( 1 )
 278 tmpdir = tempfile.mkdtemp()
 279 print("tmpdir=",tmpdir)
 280 fichierMaillage = os.path.join(tmpdir, 'garonne_1.med')
 281 garonne_1.ExportMED(fichierMaillage, 0, SMESH.MED_V2_2, 1, None ,1)
 282
 283 controlMeshStats(garonne_1, 3888, 475, 7597)
 284 controlSubMeshStats(litMineur_1, 2384)
 285 controlSubMeshStats(riveDroite_1, 2342)
 286 controlSubMeshStats(riveGauche_1, 2871, tolerance=0.1)
 287
 288 if salome.sg.hasDesktop():
 289   salome.sg.updateObjBrowser()
 290
 291 #----------------------
 292 # --- Z interpolation with HYDRO
 293 #----------------------
 294
 295 from salome.hydrotools.interpolZ import interpolZ
 296 from salome.hydrotools.controls import controlStatZ
 297
 298 # --- case name in HYDRO
 299 nomCas = 'garonne_1'
 300
 301 # --- med file 2D(x,y) of the case produced by SMESH
 302
 303 # --- dictionary [med group name] = region name
 304 dicoGroupeRegion= dict(litMineur  = 'garonne_1_litMineur',
 305                        riveDroite = 'garonne_1_riveDroite',
 306                        riveGauche = 'garonne_1_riveGauche',
 307                        )
 308 # --- value to use for Z when the node is not in a region (used to detect problems)
 309 zUndef = 90
 310 # --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation
 311 interpolMethod = 0
 312 # --- produce a 3D mesh (Z set to its value instead of 0
 313 m3d = True
 314
 315 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
 316 statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d)
 317 #print statz
 318 refstatz = {'riveDroite': (10.88, 32.61, 24.17, 5.12, 17.57, 31.33, 0.25),
 319             'riveGauche': (7.72, 71.86, 24.51, 12.18, 12.90, 60.36, 1.4),
 320             'litMineur': (2.06, 25.41, 13.93, 4.33, 8.47, 21.78)}
 321 controlStatZ(statz, refstatz)
 322
 323 # --- add a field on nodes of type double with z values, named "BOTTOM"
 324 #createZfield2(fichierMaillage)
 325 # --- create a field on nodes of type double with srtickler coefficient values, named "BOTTOM FRICTION"
 326 from salome.hydrotools.interpolS import assignStrickler
 327
 328 fichierFieldZ = os.path.join(tmpdir, 'garonne_1F.med')
 329 fichierFieldStrickler = os.path.join(tmpdir, 'garonne_1S.med')
 330
 331 assignStrickler(nomCas, fichierFieldZ, fichierFieldStrickler)