diff --git a/convert.py b/convert.py new file mode 100755 index 0000000000000000000000000000000000000000..7da196e3c2c71116dda07943bff5280b37ebae6d --- /dev/null +++ b/convert.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +# +# Quick-and-dirty script which creates a coloured .ply file from a FE calculation. +# !! works with nodal values only !! +# +# the results can be viewed with "ctmviewer" (linux) +# (sudo apt install openctm-tools) +# +# usage: +# ./fossils.py models/dolicorhynchops/dolicorhynchops_10k.py +# ./convert.py +# ctmviewer philippe.ply + +import colorsys +import vtk +import gmsh +import os.path +import numpy as np + +import fossils +fossils.setup_pythonpath() + +folder = os.path.join('workspace', \ + 'models_dolicorhynchops_dolicorhynchops_10k') +meshfile = os.path.join(folder, 'mesh.msh') +datafile = os.path.join(folder, 'smooth_stress_zz.msh') # <= values for colours + +# load mesh and result file into Gmsh +gmsh.initialize() +gmsh.merge(meshfile) +gmsh.merge(datafile) + +# retrieve view data as numpy arrays +views = gmsh.view.getTags() +datatype, tags, data, \ +time, numComponents = gmsh.view.getHomogeneousModelData(views[0], 0) + +print(f'datatype = {datatype}') +print(f'tags = {tags}') +print(f'data = {data}') +print(f'time = {time}') +print(f'numComponents = {numComponents}') + +tmp_file = 'mesh.tmp.vtk' +# write mesh to a file in VTK format +gmsh.write(tmp_file) + +# read exported mesh +reader = vtk.vtkUnstructuredGridReader() +reader.SetFileName(tmp_file) +reader.Update() +ugrid = reader.GetOutput() + +# add a colour array to the mesh +colors = vtk.vtkUnsignedCharArray() +colors.SetNumberOfComponents(3) +colors.SetName("Colors") +ugrid.GetPointData().AddArray(colors) +points = ugrid.GetPoints() +npts = points.GetNumberOfPoints() +# print(f'npts={npts}') +# print(f'len(data)={len(data)}') +vmin = np.min(data) +vmax = np.max(data) +# print(f'vmin={vmin}') +# print(f'vmax={vmax}') +rng = vmax-vmin + +for i, v in enumerate(data): + val = int((v-vmin)/rng*255) # compute "hue" value from 0 to 255 + r, g, b = colorsys.hsv_to_rgb(val/255, 1.0, 1.0) # convert to "rainbow colours" + colors.InsertTuple3(i, r*255, g*255, b*255) + +# extract mesh boundary +bfilter = vtk.vtkGeometryFilter() +bfilter.SetInputConnection(reader.GetOutputPort()) + +# export to .PLY +writer = vtk.vtkPLYWriter() +writer.SetFileName("philippe.ply") +writer.SetArrayName("Colors") +writer.SetInputConnection(bfilter.GetOutputPort()) +writer.Write() + +# clean tmp file +os.remove(tmp_file) + +print('done.')