From 274cff566b87dff63897ab08e825b198f80cb6c4 Mon Sep 17 00:00:00 2001 From: Paul RASCLE Date: Sat, 18 Jul 2020 17:12:42 +0200 Subject: [PATCH] mesh coordinates change and mesh cut functions --- src/HYDROTools/CMakeLists.txt | 2 + src/HYDROTools/changeCoords.py | 69 +++++++++++++++++++++ src/HYDROTools/cutMesh.py | 109 +++++++++++++++++++++++++++++++++ 3 files changed, 180 insertions(+) create mode 100644 src/HYDROTools/changeCoords.py create mode 100644 src/HYDROTools/cutMesh.py diff --git a/src/HYDROTools/CMakeLists.txt b/src/HYDROTools/CMakeLists.txt index 1074e2d9..b3b04c9e 100644 --- a/src/HYDROTools/CMakeLists.txt +++ b/src/HYDROTools/CMakeLists.txt @@ -25,6 +25,8 @@ SET(PYFILES interpolZ.py interpolS.py shapesGroups.py + cutMesh.py + changeCoords.py ) # --- rules --- diff --git a/src/HYDROTools/changeCoords.py b/src/HYDROTools/changeCoords.py new file mode 100644 index 00000000..8573a13b --- /dev/null +++ b/src/HYDROTools/changeCoords.py @@ -0,0 +1,69 @@ + +import numpy as np +import MEDLoader as ml +import medcoupling as mc + +from pyproj import CRS +from pyproj import Transformer + +def changeCoords(fileIn, fileOut, epsgIn=2154, epsgOut=2154, offsetXin=0, offsetYin=0, offsetXout=0, offsetYout=0): + """ + change the coordinates of a MED mesh file, with an optional offset reference in Input, and an optional change of system of coordinates (for instance Lambert II --> Lambert 93). The system of coordinates are given with EPSG (European Petroleum Survey Group) integer codes. The transformation uses pyproj. + + When there is a change of system of coordinates, input offset should be given in order to get the correct input coordinates before transformation. example, with an input file in LambertII and an origin offset of (500000, 6000000), to obtain the coordinates in Lambert93 without origin offset: changeCoords("medfileLambertII.med", "medfileLambert93.med", 27572, 2154, 500000, 6000000, 0, 0) + + When there is only an origin offset change, use the same code for epsgIn and epsgOut. + The operation is: Xout = Xin + offsetXin -offsetXout, Yout = Yin + offsetYin -offsetYout. + If the epsgIn and epsgOut codes are equal, the transformation of system of coordinates (pyproj) is not called. + + parameters: + fileIn: full path of input MED file + fileOut: full path of output MED file + epsgIn: integer code of the EPSG for the input system of coordinates (for instance, Lambert II is 27572) + epsgOut: integer code of the EPSG for the output system of coordinates (for instance, Lambert 93 is 2154) + offsetXin: X origin in the input system of coordinates (default 0) + offsetYin: Y origin in the input system of coordinates (default 0) + offsetXout: X origin in the output system of coordinates (default 0) + offsetYout: Y origin in the output system of coordinates (default 0) + return: + O if OK + """ + # TODO: do we need to transform the exception in return codes ? + if epsgIn != epsgOut: + crs_in = CRS.from_epsg(epsgIn) + crs_out = CRS.from_epsg(epsgOut) + transformer = Transformer.from_crs(crs_in, crs_out) + + meshMEDFileRead = ml.MEDFileMesh.New(fileIn) + coords = meshMEDFileRead.getCoords() + nb_comp = coords.getNumberOfComponents() + + vx = coords[:,0] + vy = coords[:,1] + avx = vx.toNumPyArray() + offsetXin + avy = vy.toNumPyArray() + offsetYin + + if nb_comp == 3: + vz = coords[:,2] + avz = vz.toNumPyArray() + + if epsgIn != epsgOut: + navx, navy = transformer.transform(avx, avy) + else: + navx = avx + navy = avy + navx = navx - offsetXout + navy = navy - offsetYout + + if nb_comp == 3: + navxy = np.stack((navx, navy, avz), axis=-1) + ncoords = mc.DataArrayDouble(navxy) + ncoords.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"]) + else: + navxy = np.stack((navx, navy), axis=-1) + ncoords = mc.DataArrayDouble(navxy) + ncoords.setInfoOnComponents(["X [m]", "Y [m]"]) + meshMEDFileRead.setCoords(ncoords) + + meshMEDFileRead.write(fileOut, 2) + return 0 diff --git a/src/HYDROTools/cutMesh.py b/src/HYDROTools/cutMesh.py new file mode 100644 index 00000000..56c18d15 --- /dev/null +++ b/src/HYDROTools/cutMesh.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python + +import os +import sys +import salome + +salome.salome_init() + +from HYDROPy import * + +def cutMesh(meshFileIn, meshFileOut = "", polyFile, offsetX=0, offsetY=0): + """ + In a given MED mesh, remove faces and edges contained in a polygon (shapefile). Save the result in a new file. + parameters: + meshFileIn: full path of input MED file. Coordinates should be without an origin offset of coordinates. + meshFileout: full path of output MED file (default="" : when "", output file is suffixed wit "_cut.med" + polyFile: a shapefile giving a single polygon for cut (should be in the same coordinate system as the mesh, without an origin offset of coordinates. + offsetX: local X origin for cut operation and output + offsetY: local Y origin for cut operation and output + """ + hydro_doc = HYDROData_Document.Document() + + # --- coordonnées locales + + hydro_doc.SetLocalCS( offsetX, offsetY ) + + # --- polygone decoupe + + a = os.path.splitext(polyFile) + polyName = os.path.basename(a[0]) + HYDROData_PolylineXY.ImportShapesFromFile(polyFile) + poly = hydro_doc.FindObjectByName(polyName + '_PolyXY_0') + + # --- zone immersible + + zoneImmersible = hydro_doc.CreateObject( KIND_IMMERSIBLE_ZONE ) + zoneImmersible.SetName( "zoneImmersible" ) + zoneImmersible.SetPolyline( poly ) + zoneImmersible.Update() + + # --- cas de calcul = une face simple sur le polygone + + casCalcul = hydro_doc.CreateObject( KIND_CALCULATION ) + casCalcul.SetName( "casCalcul" ) + casCalcul.SetAssignmentMode( HYDROData_CalculationCase.AUTOMATIC ) + casCalcul.AddGeometryObject( zoneImmersible ) + case_geom_group = zoneImmersible.GetGroup( 0 ) + casCalcul.AddGeometryGroup( case_geom_group ) + casCalcul.SetBoundaryPolyline( poly ) + casCalcul.Update() + + # --- export geom + + casCalcul_entry = casCalcul.Export() + print("casCalcul_entry", casCalcul_entry) + + import GEOM + from salome.geom import geomBuilder + geompy = geomBuilder.New() + + geomObj = salome.IDToObject(str(casCalcul_entry)) + + # --- decoupe maillage dans SMESH + + import SMESH, SALOMEDS + from salome.smesh import smeshBuilder + + # --- chargement maillage, changement de repère + + smesh = smeshBuilder.New() + ([Mesh], status) = smesh.CreateMeshesFromMED(meshFileIn) + Mesh.TranslateObject( Mesh, [ -offsetX, -offsetY, 0 ], 0 ) + + # --- groupes des faces et des edges contenues dans le polygone de découpe + + aCriteria = [] + aCriterion = smesh.GetCriterion(SMESH.FACE,SMESH.FT_BelongToGeom,SMESH.FT_Undefined, geomObj) + aCriteria.append(aCriterion) + aFilter = smesh.GetFilterFromCriteria(aCriteria) + aFilter.SetMesh(Mesh.GetMesh()) + decoupeFaces = Mesh.GroupOnFilter( SMESH.FACE, 'decoupeFaces', aFilter ) + + aCriteria = [] + aCriterion = smesh.GetCriterion(SMESH.EDGE,SMESH.FT_BelongToGeom,SMESH.FT_Undefined, geomObj) + aCriteria.append(aCriterion) + aFilter = smesh.GetFilterFromCriteria(aCriteria) + aFilter.SetMesh(Mesh.GetMesh()) + decoupeEdges = Mesh.GroupOnFilter( SMESH.EDGE, 'decoupeEdges', aFilter ) + + # --- suppression des groupes et de leurs mailles + + Mesh.RemoveGroupWithContents( decoupeFaces ) + Mesh.RemoveGroupWithContents( decoupeEdges ) + + # --- regeneration des edges de bord manquantes + + nbAdded, Mesh, _NoneGroup = Mesh.MakeBoundaryElements( SMESH.BND_1DFROM2D, '', '', 0, []) + + # --- enregistrement MED du maillage découpé + + if meshFileOut == "": + a = os.path.splitext(meshFileIn) + smesh.SetName(Mesh, os.path.basename(a[0])) + meshFileOut = a[0] + '_cut' + a[1] + + Mesh.ExportMED(meshFileOut,auto_groups=0,minor=40,overwrite=1,meshPart=None,autoDimension=1) + + if salome.sg.hasDesktop(): + salome.sg.updateObjBrowser() -- 2.39.2