save also .vtp

convert.py

@@ -21,13 +21,24 @@ import numpy as np
 import fossils
 fossils.setup_pythonpath()
 
+# - parameters
+testcase = 'models_Panthera_pardus_Panthera_pardus_RC_30'
+fieldname = 'smooth_stress_zz.msh'
+use_custom_range = True
+vmin = -30. # <= use custom range (panthera)
+vmax = 30.
+
+testcase = 'models_dolicorhynchops_dolicorhynchops_10k'
+fieldname = 'smooth_stress_zz.msh'
+use_custom_range = False
+
+# ======================================================================
 
 # folder = os.path.join('workspace', \
 #         'models_dolicorhynchops_dolicorhynchops_10k')
-folder = os.path.join('workspace', \
-        'models_Panthera_pardus_Panthera_pardus_RC_30')
+folder = os.path.join('workspace', testcase)
 meshfile = os.path.join(folder, 'mesh.msh')
-datafile = os.path.join(folder, 'smooth_stress_zz.msh') # <= values for colours
+datafile = os.path.join(folder, fieldname) # <= values for colours
 
 # load mesh and result file into Gmsh
 gmsh.initialize()
@@ -65,11 +76,9 @@ npts = points.GetNumberOfPoints()
 # print(f'npts={npts}')
 # print(f'len(data)={len(data)}')
 
-vmin = np.min(data) # <= use real full range (dolico)
-vmax = np.max(data)
-
-vmin = -30. # <= use custom range (panthera)
-vmax = 30.
+if not use_custom_range:
+    vmin = np.min(data) # <= use real full range (dolico)
+    vmax = np.max(data)
 
 # print(f'vmin={vmin}')
 # print(f'vmax={vmax}')
@@ -93,7 +102,41 @@ writer.SetArrayName("Colors")
 writer.SetInputConnection(bfilter.GetOutputPort())
 writer.Write()
 
+# ==============================================================================
+
+# ----- export .vtk with a single nodal field
+
+# read exported mesh
+reader = vtk.vtkUnstructuredGridReader()
+reader.SetFileName(tmp_file)
+reader.Update()
+ugrid = reader.GetOutput()
+
+
+darray = vtk.vtkDoubleArray()
+darray.SetNumberOfComponents(1)
+darray.SetName("sigzz")
+ugrid.GetPointData().AddArray(darray)
+ugrid.GetPointData().SetActiveScalars("sigzz")
+for i, v in enumerate(data):
+    if v>vmax: v=vmax # saturate
+    if v<vmin: v=vmin
+    darray.InsertValue(i, v)
+
+# extract mesh boundary
+bfilter = vtk.vtkGeometryFilter()
+bfilter.SetInputConnection(reader.GetOutputPort())
+
+
+writer = vtk.vtkXMLPolyDataWriter()
+writer.SetInputConnection(bfilter.GetOutputPort())
+writer.SetFileName('philippe.vtp')
+writer.Write()
+
+# ===========================================================
+
 # clean tmp file
 os.remove(tmp_file)
 
+
 print('done.')