From 31b2900d46f55330d94df371683a7aaecd13d815 Mon Sep 17 00:00:00 2001
From: mzn <mzn@opencascade.com>
Date: Fri, 23 Sep 2016 13:53:57 +0300
Subject: [PATCH] Lot5: add med field creation script.

---
 src/HYDROTools/interpolS.py | 89 +++++++++++++++++++++++++++++++++++++
 1 file changed, 89 insertions(+)
 create mode 100644 src/HYDROTools/interpolS.py

diff --git a/src/HYDROTools/interpolS.py b/src/HYDROTools/interpolS.py
new file mode 100644
index 0000000..000df60
--- /dev/null
+++ b/src/HYDROTools/interpolS.py
@@ -0,0 +1,89 @@
+# -*- coding: utf-8 -*-
+
+import os
+import sys
+import shutil
+
+import salome
+
+from med import medfile
+from med import medmesh
+from med import medfield
+from med import medenum
+
+import HYDROPy
+
+# Get current study id
+salome.salome_init()
+theStudy = salome.myStudy
+theStudyId = salome.myStudyId
+
+"""Bad parameters exception"""
+class BadParamsError(ValueError):
+  pass
+
+"""
+Creates MED file with scalar field associated with the mesh of the calculation case. 
+The scalar value associated with the mesh node is the Strickler coefficient of the land cover containing the node.
+case_name: calculation case name in the study
+med_file_name: path to input MED file with mesh on nodes
+output_file_name: path to output MED file with
+med_field_name: field name
+"""
+def assignStrickler(case_name, med_file_name, output_file_name, med_field_name='FRICTION'):
+  # Check calculation case
+  doc = HYDROPy.HYDROData_Document.Document( theStudyId )
+  case = doc.FindObjectByName(case_name)
+  if case is None:
+	raise BadParamsError("Calculation case '%s' not found" % case_name)
+	
+  # Check input MED file
+  if not os.path.exists(med_file_name):
+	raise BadParamsError("Input file '%s' not exists" % med_file_name)
+
+  # Copy input file to the output path, if the output path is empty - use the input path for the output
+  file_path = med_file_name
+  if output_file_name and (output_file_name != file_path):
+	shutil.copyfile(med_file_name, output_file_name)
+	file_path = output_file_name
+
+  # Open input MED file
+  fid = medfile.MEDfileOpen(file_path, medenum.MED_ACC_RDEXT)
+
+  # Get mesh info
+  mesh_name, sdim, mdim, meshtype, desc, dtunit, sort, nstep,  repere, axisname, axisunit = medmesh.MEDmeshInfo(fid, 1)
+  
+  nb_nodes, chgt, trsf = medmesh.MEDmeshnEntity(fid, mesh_name, medenum.MED_NO_DT, medenum.MED_NO_IT,
+												medenum.MED_NODE, medenum.MED_NONE,
+												medenum.MED_COORDINATE, medenum.MED_NO_CMODE)
+  
+  
+  # Create field
+  if nb_nodes > 0:
+    # Get node coordinates    
+    coords = medfile.MEDFLOAT(nb_nodes*sdim)
+    medmesh.MEDmeshNodeCoordinateRd(fid, mesh_name, medenum.MED_NO_DT, medenum.MED_NO_IT,
+									medenum.MED_FULL_INTERLACE, coords)
+  
+    # Get list of Strickler coefficient values for the nodes
+    values = []
+    x_coords = coords[0::sdim]
+    y_coords = coords[1::sdim]
+    values = case.GetStricklerCoefficientForPoints(x_coords, y_coords, 0.0, True)
+  
+    # Write the values to the field	
+    if len(values) > 0:
+      values = medfile.MEDFLOAT(values)
+
+      comp  = "strickler coeff "
+      medfield.MEDfieldCr(fid, med_field_name, medfile.MED_FLOAT64, 
+                          1, comp, '', '', mesh_name)
+					  
+      medfield.MEDfieldValueWr(fid, med_field_name, medenum.MED_NO_DT, medenum.MED_NO_IT, 0.0,
+                               medenum.MED_NODE, medenum.MED_NONE,
+                               medenum.MED_FULL_INTERLACE, medenum.MED_ALL_CONSTITUENT, len(values), values)
+							   
+      print "MED field '%s' on %s nodes has been created." % (med_field_name, nb_nodes)
+							 
+  # close MED file
+  medfile.MEDfileClose(fid)
\ No newline at end of file
-- 
2.39.2