4
$\begingroup$

I have a VTK file that contains several vertex fields, and I'd like to manipulate the view.

The objective is to show/hide all vertices whose field match a certain function.

For instance, if I have a temperature field, I'd like to show (all triangles) with T > -20, or even something useless like Sin[T*Norm[NormalVector]] < 0.5. However, I have no idea on how to achieve this, or how to even access one particular field in the VTK file.

Any hint is greatly appreciated!

$\endgroup$
2
  • 1
    $\begingroup$ The current version of the VTK Import, I think, only reads legacy format and ignores vertex fields. $\endgroup$
    – Tim Laska
    Dec 4 '20 at 3:00
  • $\begingroup$ @TimLaska thanks for your comment. As far as I see, in "Import Elements" it is possible to use the PointData option, with (quoting) "point data in indexed form". Maybe there is a possibility there? $\endgroup$
    – senseiwa
    Dec 5 '20 at 8:39
6
$\begingroup$

As suggested in the comments, the current version of Mathematica (12.1.1) import functionality of VTK files seems restricted to legacy format files, and it appears to ignore field data.

Sample VTK file

Here's a simple example of a cube with scaler nodal data taken from here.

SetDirectory[NotebookDirectory[]]
writeScript[filename_, txt_] := Module[{file},
  file = OpenWrite[filename];
  WriteString[file, txt];
  Close[file]]
runScript[] := RunProcess["run_script.bat"]
(*Legacy VTK format cube example*)
cubevtk = "# vtk DataFile Version 2.0
  Cube example
  ASCII
  DATASET POLYDATA
  POINTS 8 float
  0.0 0.0 0.0
  1.0 0.0 0.0
  1.0 1.0 0.0
  0.0 1.0 0.0
  0.0 0.0 1.0
  1.0 0.0 1.0
  1.0 1.0 1.0
  0.0 1.0 1.0
  POLYGONS 6 30
  4 0 1 2 3
  4 4 5 6 7
  4 0 1 5 4
  4 2 3 7 6
  4 0 4 7 3
  4 1 2 6 5
  CELL_DATA 6
  SCALARS cell_scalars int 1
  LOOKUP_TABLE default
  0
  1
  2
  3
  4
  5
  NORMALS cell_normals float
  0 0 -1
  0 0 1
  0 -1 0
  0 1 0
  -1 0 0
  1 0 0
  FIELD FieldData 2
  cellIds 1 6 int
  0 1 2 3 4 5
  faceAttributes 2 6 float
  0.0 1.0 1.0 2.0 2.0 3.0 3.0 4.0 4.0 5.0 5.0 6.0
  POINT_DATA 8
  SCALARS sample_scalars float 1
  LOOKUP_TABLE my_table
  0.0
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  LOOKUP_TABLE my_table 8
  0.0 0.0 0.0 1.0
  1.0 0.0 0.0 1.0
  0.0 1.0 0.0 1.0
  1.0 1.0 0.0 1.0
  0.0 0.0 1.0 1.0
  1.0 0.0 1.0 1.0
  0.0 1.0 1.0 1.0
  1.0 1.0 1.0 1.0";
writeScript["cube.vtk", cubevtk]
SetEnvironment["TIMSVTK_FILE" -> "cube.vtk"]

We can use a post-processor (e.g., EnSight) to view the cube.VTK file and color it by the samples_scalars variable as shown in the following image.

EnSight image

Mathematica VTK import

The following code shows which elements Mathematica imports.

elems = ImportString[cubevtk, "Elements"];
data = {Text[#], ImportString[cubevtk, #]} & /@ elems;
Text@Grid[Prepend[data, {"Element", "Returned"}], 
  Background -> {None, {Lighter[Yellow, .9], {White, 
      Lighter[Blend[{Blue, Green}], .8]}}}, 
  Dividers -> {{Darker[Gray, .6], {Lighter[Gray, .5]}, 
     Darker[Gray, .6]}, {Darker[Gray, .6], Darker[Gray, .6], {False}, 
     Darker[Gray, .6]}}, Alignment -> {{Left, Right, {Left}}}, 
  Frame -> Darker[Gray, .6], ItemStyle -> 14, 
  Spacings -> {Automatic, .8}]

Import elements

Under PointData, where sample_scalars should exist, returns an empty brace. There may be a way to get Mathematica to import the scalar field, but I have not found it. Alternatively, one could use Mathematica interoperability capability and tap into one of the freely available Python libraries to extract the field data. I will show one such approach.

Approach using external Python libraries

It is fairly straightforward to install VTK and PyVista under the Anaconda Python distribution. PyVista provides a "pythonic" interface to several types of mesh formats.

Here's an approach that I threw together that works under Windows 10. It uses a Windows batch file, run_script.bat to activate the appropriate Anaconda environment and run a Python script. I chose to use environment variables to pass in the script name and arguments versus using command-line arguments in this example. The first Python script get_arr_names.py returns a list of the available arrays in the VTK file. The second Python script, get_arr_values.py, will write an output file of the scalar values based on the array name.

winbatch = 
  Uncompress[
   "1:eJx9jzEPAUEQhfkfilFoBeVVRBQKidgTzSSyZkdu4m5Hdhfn37s7RKd975s37w1P\
ujODfq83ZyoU9HzGyXSWZXAQ7/QR4WQTFRApyDUB2bIE0qqy3oFEKOXCYBsmNojeAvHX/\
aSYt5gKBsdRAjtYeEvqnQX2dwnqK/\
appbvsZYa4jxwiYi4V4heeIZquQmNYSnK3icdNtT8H7wKrmumWft+3z1So/+\
xpmWsnjLlmGOXrjTma5W69zUetx7V0zAvEmGVI"];
writeScript["run_script.bat", winbatch];

getarrnames = 
  Uncompress[
   "1:eJxTTMoPCp7FwMCQmVuQX1SiUFBZlllckqiQWKxQUBZjYGgEFc8vBnFAuKwkOy0z\
J1XBFiiml5pXllmUnxetFOLpGxwW4h3v5unjqhQLUpebWpwBVFRQpleUmpiiAdWmCZIqLk\
3KSk0uKQZKg1TpJRYVJVbG5yUCOSDpgqLMvBINmCJNAHwpN6w="];
writeScript["get_arr_names.py", getarrnames];

getarrvals = 
  Uncompress[
   "1:eJxtj8FOwzAMhsebmJ5agaLBncMOIE3ApZuQUKmqtKRaBk0iJ43oc/\
GCOJ077cAhku3/8/87160td79Xq5UenMUAZhzcBNKDcR/ru3ueWn/\
RuClqH2SCXLyYo0pNejF89fpbwQMtCmWiRmuqbL993b3tn5un7ctjVjMnEf/\
Bmk1Zbt5PkB0DmzF+\
A5mgYbakDcofSHVRoJKfOYcXSfJje1Rd8CQnStC6nBojqUkyKtseSUMlOjs42spPGcXiTX\
WiianYC3qLsNTawDlC9zD7Ca8kdoechaKu1jVbLWc4q01o5mN8xREz45CEnH88H2Gc8DKq\
8HMe38Jy4B+cC5g4"];
writeScript["get_arr_values.py", getarrvals];

First, we will get a list of the available array names. Of which, we see sample_scalars is available.

SetEnvironment["TIMS_SCRIPT" -> "get_arr_names.py"]; 
runScript[]["StandardOutput"]
(* "['sample_scalars', 'cell_scalars', 'cell_normals', \
'cellIds', 'faceAttributes']
" *)

Now, we can write sample_scalars to an output file and read into Mathematica.

(*Write sample_scalars to output file and read into MMA*)
SetEnvironment[{"TIMS_SCRIPT" -> "get_arr_values.py", 
  "TIMS_ARRAY" -> "sample_scalars"}]
scalarfile = StringTrim@runScript[]["StandardOutput"];
scalarvals = Flatten@Import[scalarfile];
(*Read in vertex data*)
vertices = ImportString[cubevtk, "VertexData"];
{min, max} = MinMax@vertices;

We can create an interpolation function based on the VertexData and the sample_scalars. The image matches closely to what we obtained with EnSight.

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[vertices];
interp = ElementMeshInterpolation[{mesh}, scalarvals];
mr = MeshRegion@ToElementMesh[Cuboid[min {1, 1, 1}, max {1, 1, 1}]];
SliceContourPlot3D[
 interp[x, y, z], mr, {x, min, max}, {y, min, max}, {z, min, 
  max}, Contours -> 17, ColorFunction -> "ThermometerColors"]

Mathematica interpolation

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.