doc/salome/examples/h021_meshChangeBathy.py

   1 # -*- coding: utf-8 -*-
   2
   3 import os
   4 HYDRO_SAMPLES = os.path.join( os.environ["HYDRO_ROOT_DIR"], "bin/salome/test/HYDRO")
   5
   6 import sys
   7 import salome
   8
   9 salome.salome_init()
  10
  11 #----------------------
  12 # --- HYDRO
  13 #----------------------
  14
  15 from HYDROPy import *
  16 from PyQt5.QtCore import *
  17 from PyQt5.QtGui import *
  18
  19 hydro_doc = HYDROData_Document.Document()
  20
  21 hydro_doc.SetLocalCS( 430000, 6.35e+06 )
  22
  23 name = "garonne"
  24 shape = os.path.join(HYDRO_SAMPLES, name+".shp" )
  25 HYDROData_PolylineXY.ImportShapesFromFile(shape)
  26 garonne = hydro_doc.FindObjectByName(name + '_PolyXY_0')
  27 for i in range(garonne.NbSections()):
  28     garonne.SetSectionType(i, 1) # spline
  29     garonne.Update()
  30 garonne.SetZLevel( 2 )
  31
  32 name = "lit_majeur"
  33 shape = os.path.join(HYDRO_SAMPLES, name+".shp" )
  34 HYDROData_PolylineXY.ImportShapesFromFile(shape)
  35 lit_majeur = hydro_doc.FindObjectByName(name + '_PolyXY_0')
  36 for i in range(lit_majeur.NbSections()):
  37     lit_majeur.SetSectionType(i, 1) # spline
  38     lit_majeur.Update()
  39 lit_majeur.SetZLevel( 3 )
  40
  41 name = "domaine"
  42 shape = os.path.join(HYDRO_SAMPLES, name+".shp" )
  43 HYDROData_PolylineXY.ImportShapesFromFile(shape)
  44 domaine = hydro_doc.FindObjectByName(name + '_PolyXY_0')
  45 for i in range(domaine.NbSections()):
  46     domaine.SetSectionType(i, 0) # polyline
  47     domaine.Update()
  48 domaine.SetZLevel( 4 )
  49
  50 Cloud_02 = hydro_doc.CreateObject( KIND_BATHYMETRY )
  51 Cloud_02.SetName( "Cloud_02" )
  52
  53 Cloud_02.SetAltitudesInverted( 0 )
  54 if not(Cloud_02.ImportFromFile( os.path.join(HYDRO_SAMPLES, "Cloud_02.xyz" ))):
  55   raise ValueError('problem while loading bathymetry')
  56
  57 Cloud_02.Update()
  58
  59
  60 garonne_point_L93 = hydro_doc.CreateObject( KIND_BATHYMETRY )
  61 garonne_point_L93.SetName( "garonne_point_L93" )
  62
  63 garonne_point_L93.SetAltitudesInverted( 0 )
  64 if not(garonne_point_L93.ImportFromFile( os.path.join(HYDRO_SAMPLES, "garonne_point_L93.xyz" ))):
  65   raise ValueError('problem while loading bathymetry')
  66
  67 garonne_point_L93.Update()
  68
  69
  70 litMineur = hydro_doc.CreateObject( KIND_IMMERSIBLE_ZONE )
  71 litMineur.SetName( "litMineur" )
  72
  73 litMineur.SetZLevel( 6 )
  74
  75 litMineur.SetAltitudeObject( garonne_point_L93 )
  76 litMineur.SetPolyline( garonne )
  77
  78 litMineur.Update()
  79
  80
  81 litMajeur = hydro_doc.CreateObject( KIND_IMMERSIBLE_ZONE )
  82 litMajeur.SetName( "litMajeur" )
  83
  84 litMajeur.SetZLevel( 5 )
  85
  86 litMajeur.SetFillingColor( QColor( 0, 170, 127, 255 ) )
  87
  88 litMajeur.SetAltitudeObject( garonne_point_L93 )
  89 litMajeur.SetPolyline( lit_majeur )
  90
  91 litMajeur.Update()
  92
  93
  94 domaineEtendu = hydro_doc.CreateObject( KIND_IMMERSIBLE_ZONE )
  95 domaineEtendu.SetName( "domaineEtendu" )
  96
  97 domaineEtendu.SetZLevel( 4 )
  98
  99 domaineEtendu.SetFillingColor( QColor( 201, 203, 55, 255 ) )
 100
 101 domaineEtendu.SetAltitudeObject( Cloud_02 )
 102 domaineEtendu.SetPolyline( domaine )
 103
 104 domaineEtendu.Update()
 105
 106
 107 # Calculation case
 108 garonne_1 = hydro_doc.CreateObject( KIND_CALCULATION )
 109 garonne_1.SetName( "garonne_1" )
 110
 111 garonne_1.SetAssignmentMode( HYDROData_CalculationCase.MANUAL )
 112 garonne_1.AddGeometryObject( litMineur )
 113 garonne_1.AddGeometryObject( domaineEtendu )
 114 garonne_1.AddGeometryObject( litMajeur )
 115
 116 case_geom_group = domaineEtendu.GetGroup( 0 )
 117 garonne_1.AddGeometryGroup( case_geom_group )
 118 case_geom_group = litMineur.GetGroup( 0 )
 119 garonne_1.AddGeometryGroup( case_geom_group )
 120 case_geom_group = litMajeur.GetGroup( 0 )
 121 garonne_1.AddGeometryGroup( case_geom_group )
 122 garonne_1.SetBoundaryPolyline( domaine )
 123
 124 # Start the algorithm of the partition and assignment
 125 garonne_1.Update()
 126 garonne_1_litMineur = hydro_doc.FindObjectByName( "garonne_1_Reg_1" )
 127 garonne_1_Zone_1 = hydro_doc.FindObjectByName( "garonne_1_Zone_1" )
 128 garonne_1_Zone_1.SetMergeType( HYDROData_Zone.Merge_ZMIN )
 129 garonne_1_Zone_1.SetColor( QColor( 192, 113, 64 ))
 130 garonne_1_litMineur.AddZone( garonne_1_Zone_1 )
 131
 132 garonne_1_riveDroite = hydro_doc.FindObjectByName( "garonne_1_Reg_2" )
 133 garonne_1_Zone_2 = hydro_doc.FindObjectByName( "garonne_1_Zone_2" )
 134 garonne_1_Zone_2.SetColor( QColor( 141, 192, 64 ))
 135 garonne_1_riveDroite.AddZone( garonne_1_Zone_2 )
 136
 137 garonne_1_Zone_3 = hydro_doc.FindObjectByName( "garonne_1_Zone_3" )
 138 garonne_1_Zone_3.SetMergeType( HYDROData_Zone.Merge_Object )
 139 Cloud_02 = hydro_doc.FindObjectByName( "Cloud_02" )
 140 garonne_1_Zone_3.SetMergeObject( Cloud_02 )
 141 garonne_1_Zone_3.SetColor( QColor( 64, 192, 77 ))
 142 garonne_1_riveDroite.AddZone( garonne_1_Zone_3 )
 143
 144 garonne_1_riveGauche = hydro_doc.FindObjectByName( "garonne_1_Reg_3" )
 145 garonne_1_Zone_4 = hydro_doc.FindObjectByName( "garonne_1_Zone_4" )
 146 garonne_1_Zone_4.SetMergeType( HYDROData_Zone.Merge_Object )
 147 Cloud_02 = hydro_doc.FindObjectByName( "Cloud_02" )
 148 garonne_1_Zone_4.SetMergeObject( Cloud_02 )
 149 garonne_1_Zone_4.SetColor( QColor( 64, 75, 192 ))
 150 garonne_1_riveGauche.AddZone( garonne_1_Zone_4 )
 151
 152 garonne_1_Zone_5 = hydro_doc.FindObjectByName( "garonne_1_Zone_5" )
 153 garonne_1_Zone_5.SetColor( QColor( 64, 192, 77 ))
 154 garonne_1_riveGauche.AddZone( garonne_1_Zone_5 )
 155
 156 garonne_1_litMineur.SetName("garonne_1_litMineur")
 157 garonne_1_riveDroite.SetName("garonne_1_riveDroite")
 158 garonne_1_riveGauche.SetName("garonne_1_riveGauche")
 159
 160 # Export of the calculation case
 161 garonne_1_entry = garonne_1.Export()
 162
 163 # --- add a new bathymetry for the test changeBathy
 164
 165 import tempfile
 166 tmpdir = tempfile.mkdtemp()
 167 print("tmpdir=",tmpdir)
 168
 169 newBathy = os.path.join(tmpdir, 'newBathy.xyz')
 170 fi=open(os.path.join(HYDRO_SAMPLES, "garonne_point_L93.xyz" ), 'r')
 171 fo=open(newBathy, 'w')
 172 for ligne in fi:
 173     vals = ligne.split()
 174     if len(vals) < 3:
 175         continue
 176     x = float(vals[0])
 177     y = float(vals[1])
 178     z = float(vals[2]) + 50
 179     l = "%12.3f  %12.3f %12.3f\n" % (x, y, z)
 180     fo.write(l)
 181 fi.close()
 182 fo.close()
 183
 184 #----------------------
 185 # --- Geometry
 186 #----------------------
 187
 188 # Get geometry shape and print debug information
 189 import GEOM
 190 from salome.geom import geomBuilder
 191 import math
 192 import SALOMEDS
 193 from salome.hydrotools.controls import controlGeomProps
 194
 195 geompy = geomBuilder.New()
 196
 197 print("Entry:", garonne_1_entry)
 198 HYDRO_garonne_1 = salome.IDToObject( str( garonne_1_entry ) )
 199 print("Geom shape:", HYDRO_garonne_1)
 200 print("Geom shape name:", HYDRO_garonne_1.GetName())
 201
 202 # --- manual definition: geometrical faces
 203
 204 [garonne_litMineur,garonne_riveDroite,garonne_riveGauche] = geompy.SubShapeAll(HYDRO_garonne_1, geompy.ShapeType["FACE"])
 205
 206 controlGeomProps(geompy, garonne_riveGauche,  29149.36918,  35948828.352061)
 207 controlGeomProps(geompy, garonne_litMineur,   30337.548492,  3488480.304388)
 208 controlGeomProps(geompy, garonne_riveDroite,  32012.343241, 25998769.23615)
 209
 210 # --- manual identification of all useful edge groups (boundary conditions)
 211
 212 allEdgesIds = geompy.SubShapeAllIDs(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 213 print("allEdgesIds", allEdgesIds)
 214
 215 (isDone, ClosedFreeBoundary, OpenFreeBoundary) = geompy.GetFreeBoundary(HYDRO_garonne_1)
 216 geompy.addToStudyInFather(HYDRO_garonne_1, ClosedFreeBoundary[0], "ClosedFreeBoundary")
 217
 218 freeBoundary = geompy.ExtractShapes(ClosedFreeBoundary[0], geompy.ShapeType["EDGE"], True)
 219 freeBoundaryIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, freeBoundary[i]) for i in range(len(freeBoundary)) ]
 220 print("freeBoundaryIds", freeBoundaryIds)
 221
 222 [litMineur_droite] = geompy.GetSharedShapesMulti([garonne_riveDroite, garonne_litMineur], geompy.ShapeType["EDGE"], True)
 223 [litMineur_gauche] = geompy.GetSharedShapesMulti([garonne_riveGauche, garonne_litMineur], geompy.ShapeType["EDGE"], True)
 224 geompy.addToStudyInFather(HYDRO_garonne_1, litMineur_droite, "litMineur_droite")
 225 geompy.addToStudyInFather(HYDRO_garonne_1, litMineur_gauche, "litMineur_gauche")
 226 rives = [litMineur_droite, litMineur_gauche]
 227 rivesIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, rives[i]) for i in range(len(rives)) ]
 228 print("rivesIds", rivesIds)
 229
 230 edges_litMineur = geompy.GetSharedShapesMulti([HYDRO_garonne_1, garonne_litMineur], geompy.ShapeType["EDGE"], True)
 231 edges_riveGauche = geompy.GetSharedShapesMulti([HYDRO_garonne_1, garonne_riveGauche], geompy.ShapeType["EDGE"], True)
 232 edges_riveDroite = geompy.GetSharedShapesMulti([HYDRO_garonne_1, garonne_riveDroite], geompy.ShapeType["EDGE"], True)
 233 edges_litMineurIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_litMineur[i]) for i in range(len(edges_litMineur)) ]
 234 edges_riveGaucheIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveGauche[i]) for i in range(len(edges_riveGauche)) ]
 235 edges_riveDroiteIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveDroite[i]) for i in range(len(edges_riveDroite)) ]
 236
 237 print("edges_litMineurIds", edges_litMineurIds)
 238 print("edges_riveGaucheIds", edges_riveGaucheIds)
 239 print("edges_riveDroiteIds", edges_riveDroiteIds)
 240
 241 sectionsIds = [Id for Id in edges_litMineurIds if Id not in rivesIds]
 242 print("sectionsIds", sectionsIds)
 243 SectionsGaronne = geompy.CreateGroup(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 244 geompy.UnionIDs(SectionsGaronne, sectionsIds)
 245 geompy.addToStudyInFather(HYDRO_garonne_1, SectionsGaronne, "SectionsGaronne")
 246
 247 bordGaucheDomaineIds = [Id for Id in freeBoundaryIds if Id in edges_riveGaucheIds]
 248 bordDroiteDomaineIds = [Id for Id in freeBoundaryIds if Id in edges_riveDroiteIds]
 249 print("bordGaucheDomaineIds", bordGaucheDomaineIds)
 250 print("bordDroiteDomaineIds", bordDroiteDomaineIds)
 251 bordGaucheDomaine = geompy.CreateGroup(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 252 geompy.UnionIDs(bordGaucheDomaine, bordGaucheDomaineIds)
 253 geompy.addToStudyInFather(HYDRO_garonne_1, bordGaucheDomaine, "bordGaucheDomaine")
 254 bordDroiteDomaine = geompy.CreateGroup(HYDRO_garonne_1, geompy.ShapeType["EDGE"])
 255 geompy.UnionIDs(bordDroiteDomaine, bordDroiteDomaineIds)
 256 geompy.addToStudyInFather(HYDRO_garonne_1, bordDroiteDomaine, "bordDroiteDomaine")
 257
 258 amont = geompy.GetEdgeNearPoint(HYDRO_garonne_1, geompy.MakeVertex(46757.861314, 25833.234752, 0))
 259 aval = geompy.GetEdgeNearPoint(HYDRO_garonne_1, geompy.MakeVertex(39078.979127, 32588.627279, 0))
 260 geompy.addToStudyInFather(HYDRO_garonne_1, amont, "amont")
 261 geompy.addToStudyInFather(HYDRO_garonne_1, aval, "aval")
 262
 263 #----------------------
 264 # --- Meshing
 265 #----------------------
 266
 267 import  SMESH, SALOMEDS
 268 from salome.smesh import smeshBuilder
 269 from salome.hydrotools.controls import controlMeshStats, controlSubMeshStats
 270 import tempfile
 271
 272 smesh = smeshBuilder.New()
 273
 274 # --- algorithms and hypothesis
 275 garonne_1 = smesh.Mesh(HYDRO_garonne_1)
 276
 277 NETGEN_2D = garonne_1.Triangle(algo=smeshBuilder.NETGEN_1D2D)
 278 NETGEN_2D_Parameters = NETGEN_2D.Parameters()
 279 NETGEN_2D_Parameters.SetMaxSize( 200 )
 280 NETGEN_2D_Parameters.SetSecondOrder( 0 )
 281 NETGEN_2D_Parameters.SetOptimize( 1 )
 282 NETGEN_2D_Parameters.SetFineness( 4 )
 283 NETGEN_2D_Parameters.SetMinSize( 50 )
 284 NETGEN_2D_Parameters.SetUseSurfaceCurvature( 1 )
 285 NETGEN_2D_Parameters.SetFuseEdges( 1 )
 286 NETGEN_2D_Parameters.SetQuadAllowed( 0 )
 287
 288 algo2D_litMineur = garonne_1.Quadrangle(algo=smeshBuilder.QUAD_MA_PROJ,geom=garonne_litMineur)
 289 algo1D_litMineur = garonne_1.Segment(geom=garonne_litMineur)
 290 hypo1D_litMineur = algo1D_litMineur.LocalLength(100,None,1e-07)
 291 subMesh_litMineur = algo1D_litMineur.GetSubMesh()
 292 smesh.SetName(subMesh_litMineur, "litMineur")
 293
 294 algo1D_SectionsGaronne = garonne_1.Segment(geom=SectionsGaronne)
 295 hypo1D_SectionsGaronne = algo1D_SectionsGaronne.NumberOfSegments(8)
 296 hypo1D_SectionsGaronne.SetDistrType( 0 )
 297 subMesh_SectionsGaronne = algo1D_SectionsGaronne.GetSubMesh()
 298 smesh.SetName(subMesh_SectionsGaronne, "SectionsGaronne")
 299
 300 isDone = garonne_1.SetMeshOrder( [ [ subMesh_SectionsGaronne, subMesh_litMineur ] ])
 301
 302 # --- compute mesh
 303 isDone = garonne_1.Compute()
 304 isDone = garonne_1.SplitQuadObject( garonne_1, 1 )
 305 isDone = garonne_1.ReorientObject( garonne_1 )
 306
 307 # --- geometrical groups of faces
 308 riveGauche_1 = garonne_1.GroupOnGeom(garonne_riveGauche,'riveGauche',SMESH.FACE)
 309 litMineur_1 = garonne_1.GroupOnGeom(garonne_litMineur,'litMineur',SMESH.FACE)
 310 riveDroite_1 = garonne_1.GroupOnGeom(garonne_riveDroite,'riveDroite',SMESH.FACE)
 311
 312 # --- geometrical groups of edges
 313
 314 ClosedFreeBoundary_1 = garonne_1.GroupOnGeom(ClosedFreeBoundary[0],'ClosedFreeBoundary',SMESH.EDGE)
 315 litMineur_droite_1 = garonne_1.GroupOnGeom(litMineur_droite,'litMineur_droite',SMESH.EDGE)
 316 litMineur_gauche_1 = garonne_1.GroupOnGeom(litMineur_gauche,'litMineur_gauche',SMESH.EDGE)
 317 SectionsGaronne_1 = garonne_1.GroupOnGeom(SectionsGaronne,'SectionsGaronne',SMESH.EDGE)
 318 bordGaucheDomaine_1 = garonne_1.GroupOnGeom(bordGaucheDomaine,'bordGaucheDomaine',SMESH.EDGE)
 319 bordDroiteDomaine_1 = garonne_1.GroupOnGeom(bordDroiteDomaine,'bordDroiteDomaine',SMESH.EDGE)
 320 amont_1 = garonne_1.GroupOnGeom(amont,'amont',SMESH.EDGE)
 321 aval_1 = garonne_1.GroupOnGeom(aval,'aval',SMESH.EDGE)
 322
 323 # --- geometrical groups of nodes
 324
 325 garonne_1_litMineur_2 = garonne_1.GroupOnGeom(garonne_litMineur,'garonne_1_litMineur',SMESH.NODE)
 326 garonne_1_riveDroite_2 = garonne_1.GroupOnGeom(garonne_riveDroite,'garonne_1_riveDroite',SMESH.NODE)
 327 garonne_1_riveGauche_2 = garonne_1.GroupOnGeom(garonne_riveGauche,'garonne_1_riveGauche',SMESH.NODE)
 328 ClosedFreeBoundary_2 = garonne_1.GroupOnGeom(ClosedFreeBoundary[0],'ClosedFreeBoundary',SMESH.NODE)
 329 litMineur_droite_2 = garonne_1.GroupOnGeom(litMineur_droite,'litMineur_droite',SMESH.NODE)
 330 litMineur_gauche_2 = garonne_1.GroupOnGeom(litMineur_gauche,'litMineur_gauche',SMESH.NODE)
 331 SectionsGaronne_2 = garonne_1.GroupOnGeom(SectionsGaronne,'SectionsGaronne',SMESH.NODE)
 332 bordGaucheDomaine_2 = garonne_1.GroupOnGeom(bordGaucheDomaine,'bordGaucheDomaine',SMESH.NODE)
 333 bordDroiteDomaine_2 = garonne_1.GroupOnGeom(bordDroiteDomaine,'bordDroiteDomaine',SMESH.NODE)
 334 amont_2 = garonne_1.GroupOnGeom(amont,'amont',SMESH.NODE)
 335 aval_2 = garonne_1.GroupOnGeom(aval,'aval',SMESH.NODE)
 336
 337 garonne_1.SetAutoColor( 1 )
 338 tmpdir = tempfile.mkdtemp()
 339 print("tmpdir=",tmpdir)
 340 fichierMaillage = os.path.join(tmpdir, 'garonne_1.med')
 341 garonne_1.ExportMED(fichierMaillage, 0, SMESH.MED_V2_2, 1, None ,1)
 342
 343 controlMeshStats(garonne_1, 3888, 475, 7597)
 344 controlSubMeshStats(litMineur_1, 2384)
 345 controlSubMeshStats(riveDroite_1, 2342)
 346 controlSubMeshStats(riveGauche_1, 2871)
 347
 348 if salome.sg.hasDesktop():
 349   salome.sg.updateObjBrowser()
 350
 351 #----------------------
 352 # --- Z interpolation with HYDRO
 353 #----------------------
 354
 355 from salome.hydrotools.interpolZ import interpolZ
 356 from salome.hydrotools.controls import controlStatZ
 357
 358 # --- case name in HYDRO
 359 nomCas = 'garonne_1'
 360
 361 # --- med file 2D(x,y) of the case produced by SMESH
 362
 363 # --- dictionary [med group name] = region name
 364 dicoGroupeRegion= dict(litMineur  = 'garonne_1_litMineur',
 365                        riveDroite = 'garonne_1_riveDroite',
 366                        riveGauche = 'garonne_1_riveGauche',
 367                        )
 368 # --- value to use for Z when the node is not in a region (used to detect problems)
 369 zUndef = 90
 370 # --- interpolation Method: 0 = nearest point on bathymetry (default), 1 = linear interpolation
 371 interpolMethod = 0
 372 # --- produce a 3D mesh (Z set to its value instead of 0
 373 m3d = True
 374
 375 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
 376 statz = interpolZ(nomCas, fichierMaillage, dicoGroupeRegion, zUndef, interpolMethod, m3d)
 377 #print statz
 378 refstatz = {'riveDroite': (10.88, 32.61, 24.17, 5.12, 17.57, 31.33, 0.25),
 379             'riveGauche': (7.72, 71.86, 24.51, 12.18, 12.90, 60.36, 0.25),
 380             'litMineur': (2.06, 25.41, 13.93, 4.33, 8.47, 21.78)}
 381 controlStatZ(statz, refstatz)
 382
 383 # --- interpolation on a new bathymetry for 'riveGauche', without using HYDRO (no Calculation Case)
 384
 385 meshFileIn = os.path.splitext(fichierMaillage)[0] + 'F.med'
 386 meshFileOut = os.path.splitext(meshFileIn)[0] + 'F.med'
 387
 388 # --- Z interpolation on the bathymety/altimetry on the mesh nodes
 389
 390 from salome.hydrotools.changeBathy import changeBathy
 391 statz = changeBathy(meshFileIn, meshFileOut, newBathy, 'riveGauche', 430000, 6350000, stats=True)
 392 print(statz)
 393 refstatz = {'riveGauche': (57.72, 77.14, 71.33, 3.36, 62.92, 75.0)}
 394 controlStatZ(statz, refstatz)
 395