From 0f4e70091d02ecd65a23e391585d0c7b31f8ad94 Mon Sep 17 00:00:00 2001 From: Paul RASCLE Date: Mon, 24 Aug 2020 18:26:08 +0200 Subject: [PATCH] test case h021_meshChangeBathy working --- doc/salome/examples/h021_meshChangeBathy.py | 36 +++- src/HYDROTools/CMakeLists.txt | 1 + src/HYDROTools/changeBathy.py | 203 ++++++++++++++++++++ 3 files changed, 237 insertions(+), 3 deletions(-) create mode 100644 src/HYDROTools/changeBathy.py diff --git a/doc/salome/examples/h021_meshChangeBathy.py b/doc/salome/examples/h021_meshChangeBathy.py index e9213d21..b18a67e7 100644 --- a/doc/salome/examples/h021_meshChangeBathy.py +++ b/doc/salome/examples/h021_meshChangeBathy.py @@ -160,6 +160,27 @@ garonne_1_riveGauche.SetName("garonne_1_riveGauche") # Export of the calculation case garonne_1_entry = garonne_1.Export() +# --- add a new bathymetry for the test changeBathy + +import tempfile +tmpdir = tempfile.mkdtemp() +print("tmpdir=",tmpdir) + +newBathy = os.path.join(tmpdir, 'newBathy.xyz') +fi=open(os.path.join(HYDRO_SAMPLES, "garonne_point_L93.xyz" ), 'r') +fo=open(newBathy, 'w') +for ligne in fi: + vals = ligne.split() + if len(vals) < 3: + continue + x = float(vals[0]) + y = float(vals[1]) + z = float(vals[2]) + 50 + l = "%12.3f %12.3f %12.3f\n" % (x, y, z) + fo.write(l) +fi.close() +fo.close() + #---------------------- # --- Geometry #---------------------- @@ -213,7 +234,7 @@ edges_litMineurIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_litMineur[i]) edges_riveGaucheIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveGauche[i]) for i in range(len(edges_riveGauche)) ] edges_riveDroiteIds = [ geompy.GetSubShapeID(HYDRO_garonne_1, edges_riveDroite[i]) for i in range(len(edges_riveDroite)) ] -print("edges_litMineurIds", edges_litMineurIds) +print("edges_litMineurIds", edges_litMineurIds) print("edges_riveGaucheIds", edges_riveGaucheIds) print("edges_riveDroiteIds", edges_riveDroiteIds) @@ -359,7 +380,16 @@ refstatz = {'riveDroite': (10.88, 32.61, 24.17, 5.12, 17.57, 31.33, 0.25), 'litMineur': (2.06, 25.41, 13.93, 4.33, 8.47, 21.78)} controlStatZ(statz, refstatz) -# --- add a field on nodes of type double with z values, named "BOTTOM" -#createZfield2(fichierMaillage) +# --- interpolation on a new bathymetry for 'riveGauche', without using HYDRO (no Calculation Case) + +meshFileIn = os.path.splitext(fichierMaillage)[0] + 'F.med' +meshFileOut = os.path.splitext(meshFileIn)[0] + 'F.med' +# --- Z interpolation on the bathymety/altimetry on the mesh nodes + +from salome.hydrotools.changeBathy import changeBathy +statz = changeBathy(meshFileIn, meshFileOut, newBathy, 'riveGauche', 430000, 6350000, stats=True) +print(statz) +refstatz = {'riveGauche': (57.72, 77.14, 71.33, 3.36, 62.92, 75.0)} +controlStatZ(statz, refstatz) diff --git a/src/HYDROTools/CMakeLists.txt b/src/HYDROTools/CMakeLists.txt index 013d25b6..49232045 100644 --- a/src/HYDROTools/CMakeLists.txt +++ b/src/HYDROTools/CMakeLists.txt @@ -27,6 +27,7 @@ SET(PYFILES shapesGroups.py cutMesh.py changeCoords.py + changeBathy.py ) # --- plugins dialogs diff --git a/src/HYDROTools/changeBathy.py b/src/HYDROTools/changeBathy.py new file mode 100644 index 00000000..d9285d69 --- /dev/null +++ b/src/HYDROTools/changeBathy.py @@ -0,0 +1,203 @@ +# -*- coding: utf-8 -*- + +# ----------------------------------------------------------------------------- + +import numpy as np +import MEDLoader as ml +from scipy import spatial + + + +def changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX=0., offsetY=0., zUndef=-100., interpMethod=0, m3d=False, stats=False): + """ + ChangeBathy takes a 2D (x,y) mesh with an already existing field of double ('BOTTOM') for altimetry. + The mesh may already have a z coordinate set with the z altimetry. + The new altimetry from the bathymetry file is set on the nodes of the selected group of face. + + In: + meshFileIn: med file name produced by SMESH, corresponding to the HYDRO case + bathyFile: bathymetry file + groupName: name of the group of faces to modify. + offsetX, offsetY: coordinates of origin in MESH + zUndef: Z value to use for nodes outside the bathymetry interpolation (should not happen: if so, algoritm problem). + default value is -100. + interpMethod: integer value + 0 = nearest point on bathymetry + 1 = linear interpolation + m3d: True/False to produce a 3D mesh. Default is False. + stats: compute statistics for SALOME test + Out: + statz: statistique for z + Key: face group name + Value: (minz, maxz, meanz, stdz, v05z, v95z) + Out: + meshFileOut : med file with Z value in a field "BOTTOM" + Option: Z value also in Z coordinate if m3D is true + """ + + # --- load the bathymetry in a KD tree, working with the same origin coordinates as the mesh + + fi=open(bathyFile, 'r') + + minx = 1.e30 + maxx = -1.e30 + miny = 1.e30 + maxy = -1.e30 + minz = 1.e30 + maxz = -1.e30 + + coords = [] + alti = [] + for ligne in fi: + vals = ligne.split() + if len(vals) < 3: + continue + x = float(vals[0]) - offsetX + y = float(vals[1]) - offsetY + z = float(vals[2]) + coords.append((x, y)) + alti.append(z) + minx = min(x, minx) + miny = min(y, miny) + minz = min(z, minz) + maxx = max(x, maxx) + maxy = max(y, maxy) + maxz = max(z, maxz) + fi.close() + + print("bornes x", minx, maxx) + print("bornes y", miny, maxy) + print("bornes z", minz, maxz) + print("nombre points:", len(alti)) + + bathyPts = np.array(coords) + print(bathyPts.shape) + + tree = spatial.KDTree(bathyPts) + + # --- load the nodes coordinates of the mesh group + + mcMesh = ml.MEDFileMesh.New(meshFileIn) + dim = mcMesh.getMeshDimension() + d0 = 0 # when dimension 2, faces are dim + if dim == 3: # when dimension 3, edges are dim -1 + d0=-1 + + groupNames = mcMesh.getGroupsOnSpecifiedLev(d0) #names of face groups + + if groupName not in groupNames: + print("there is no group of faces named %s in the mesh" %groupName) + #return + + group = mcMesh.getGroup(d0, groupName) + + # nodes ids in use: number of nodes and order of the mesh, nodes not in the group marked -1, + # nodes indices = 0 to nbNodeInGroup, not their indice in mesh! + + AllNodeIds = group.getNodeIdsInUse()[0].toNumPyArray() + print("AllNodeIds.size:", AllNodeIds.size) + cond = AllNodeIds != -1 + #nodeIds = np.extract(cond, AllNodeIds) + nids = [i for i in range(AllNodeIds.size) if int(AllNodeIds[i]) >=0] # indices in mesh + nodeIds = np.array(nids) + print("nodeIds.size:", nodeIds.size) + print("nodeIds:", nodeIds) + + allCoords = group.getCoords().toNumPyArray() + print("allCoords.shape:",allCoords.shape) + dimCoord = allCoords.shape[1] + if dimCoord == 3: + condition = np.vstack((cond, cond, cond)).transpose() + else: + condition = np.vstack((cond, cond)).transpose() + coords = np.extract(condition, allCoords) + nbRaws = coords.size // dimCoord + coords = coords.reshape(nbRaws, dimCoord) + print("coords.shape:",coords.shape) + + if coords.shape[1] == 3: + coordsXY = coords[:, 0:2] # 2 first columns + coordsZ = coords[:, 2] # 3rd column + else: + coordsXY = coords + coordsZ = None + + if interpMethod != 0: + print("TODO: the linear interpolation on bathymetry is not yet implemented, using nearest point") + + # --- get nearest node Ids for nodes of group + + (distances, ids) = tree.query(coordsXY) + print("ids.size", ids.size) + print("ids", ids) + + # --- check if there is a 'BOTTOM' Field + + try: + bottomField = ml.ReadField(meshFileIn, 'BOTTOM') + bottomArray = bottomField.getArray() + except ml.InterpKernelException: + print("no BOTTOM array") + bottomField = None + + if bottomField is None: + bottomField = ml.MEDCouplingFieldDouble(ml.ON_NODES) + bottomField.setName('BOTTOM') + bottomField.setMesh(mcMesh.getMeshAtLevel(0)) + npField = np.zeros((mcMesh.getNumberOfNodes())) + bottomArray = ml.DataArrayDouble(npField) + + # --- set altitude on field + + vz = [] + for i,id in enumerate(ids): + #print(i, "---", nodeIds[i], "---", id, alti[id]) + j = int(nodeIds[i]) + bottomArray[j] = alti[id] + if stats: + vz.append(alti[id]) + + bottomField.setArray(bottomArray) + bottomField.setTime(0.0, 0, -1) + + if (coordsZ is not None) or m3d: + print("TODO: coordsZ") + mlCoords = group.getCoords() + coords3D = ml.DataArrayDouble.Meld([mlCoords[:,0:2], bottomArray]) + coords3D.setInfoOnComponents(["X [m]", "Y [m]", "Z [m]"]) + mcMesh.setCoords(coords3D) + + mcMesh.write(meshFileOut, 2) + ml.WriteFieldUsingAlreadyWrittenMesh(meshFileOut, bottomField) + + if stats: + statz = {} + minz = np.amin(vz) + maxz = np.amax(vz) + meanz = np.mean(vz) + stdz = np.std(vz) + v05z = np.percentile(vz, 0o5) + v95z = np.percentile(vz, 95) + ligne = ".. Minimum: %f" % minz + ligne += ", maximum: %f" % maxz + ligne += ", mean: %f\n" % meanz + ligne += ".. stdeviation: %f" % stdz + ligne += ", v05z: %f" % v05z + ligne += ", v95z: %f" % v95z + print (ligne) + statz[groupName] = (minz, maxz, meanz, stdz, v05z, v95z) + return statz + + +if __name__=='__main__': + meshFileIn = 'garonne_1F.med' + meshFileOut = 'garonne_1FF.med' + bathyFile = 'newBathy.xyz' + #groupName = 'litMineur' + groupName = 'riveGauche' + #groupName = 'riveDroite' + offsetX = 430000 + offsetY = 6350000 + stats = changeBathy(meshFileIn, meshFileOut, bathyFile, groupName, offsetX, offsetY, stats=True) + + -- 2.39.2