Intro to meshio
In the previous episodes we manipulated some of the result fields from Elmer with meshio, a very nice Python package that allows us to do many useful operations with meshes and fields. However, in the episodes we did not focus into the details of how to use the package. Whilst this webiste is perhaps not the best place to do so, as the best place to understand how to use a package is, of course, its official documentation, it is useful to collect here few examples and remarks about using it with Julia, the programming language we are using for this project. This will also help clarify what was done in the previous episodes.
Installing meshio
Assuming that you have up to date Julia and Python versions installed on your system, simply follow the installation instructions in the meshio documentation.
The next step is to install PyCall. Follow the instructions provided in the documentation. Thanks to PyCall we will be able to use python libraries, like meshio, in Julia.
A Different Installation Approach
Whilst the suggestion above should get most users started, what I prefer to do is to create a python virtual environment for Julia only, prep it up with all the package I need, and make sure that PyCall is built to use this environment. Here how to do so.
Create the python virtual environment and install the required packages in it, in this case meshio:
# Setting up a Python venv
python3 -m venv juliaenv
# Installing Plotly in the venv
juliaenv/bin/pip3 install meshio
The virtual environment can be in any directory and have any name. In case your change yours you will need to adapt the following instructions.
After the virtual environment is ready, open Julia and install PyCall. First, tell Julia to search for python in our newly created python environment:
julia> ENV["PYTHON"]="~/juliaenv/bin/python3"
Then, switch to package mode by pressing ] and install PyCall and build:
pkg> add PyCall
pkg> build PyCall
You can use CTRL+C to exit package mode after the operations above have concluded.
At this point, meshio will be available in your Julia installation through PyCall. To make more python packages available, simply install them in your virtual environment.
Using meshio
The best way to understand how to use meshio from Julia is perhaps with a REPL example. So let’s just open Julia and prepare our environment.
using PyCall
After we do that we are ready to import the meshio python module into Julia:
meshio = pyimport("meshio")
The REPL should output something like this:
PyObject <module 'meshio' from '/usr/lib/python3.8/site-packages/meshio/__init__.py'>
We can now use meshio to read VTU files and manipulate them. In the examples below we will use the solution field from the Rigid Walled Room Revisited - Part 2 episode, which for convenience you can download here. It is assumed that the file is saved in the directory where you launched Julia.
Loading a VTU File
Loading of files is very easy. Just do:
data = meshio.read("case_t0001.vtu")
The REPL should output something like this:
PyObject <meshio mesh object>
Number of points: 17949
Number of cells:
tetra: 94353
line: 284
triangle: 7474
Point data: pressure EigenMode1, pressure EigenMode2, pressure EigenMode3, pressure EigenMode4, pressure EigenMode5, pressure EigenMode6, pressure EigenMode7, pressure EigenMode8, pressure EigenMode9, pressure EigenMode10
Cell data: GeometryIds
The Data Structure
You will probably recognise many of the terms above as they are also displayed by ParaView (see the Intro to ParaView episode for more info). Our solver exported a number of scalar fields. Those are available as Point data.
To interact with the data object we will need to access its attributes. This is simply done with a . followed by the attributed name. In the REPL, you can type data. and then press the tab key two times to get a bunch of suggestions that should look like this:
__class__ __gt__ __reduce__ cell_data_dict info
__delattr__ __hash__ __reduce_ex__ cell_sets int_data_to_sets
__dict__ __init__ __repr__ cell_sets_dict point_data
__dir__ __init_subclass__ __setattr__ cells point_sets
__doc__ __le__ __sizeof__ cells_dict points
__eq__ __lt__ __str__ field_data prune
__format__ __module__ __subclasshook__ get_cell_data read
__ge__ __ne__ __weakref__ get_cells_type sets_to_int_data
__getattribute__ __new__ cell_data gmsh_periodic write
This lists all attributes and methods of the data object. We are mainly interested in point_data, as our eigenmode scalar fields are all there. So let’s take it:
field_data = data.point_data
As you can see from the REPL, the point_data attribute is a dictionary:
Dict{Any,Any} with 10 entries:
"pressure EigenMode5" => [-0.988624; 0.986744; … ; -0.0672767; 0.280671]
"pressure EigenMode9" => [-1.0; 0.996853; … ; -0.00439771; 0.0217661]
"pressure EigenMode3" => [-0.996737; -0.996578; … ; 0.446996; 0.329133]
"pressure EigenMode2" => [0.998803; 0.997648; … ; -0.138979; -0.234202]
"pressure EigenMode8" => [0.995937; -0.994356; … ; -0.0314515; 0.0934631]
"pressure EigenMode1" => [-1.0; -1.0; … ; -1.0; -1.0]
"pressure EigenMode10" => [-0.994816; -0.991452; … ; -0.427958; -0.291705]
"pressure EigenMode6" => [0.997145; -0.994048; … ; -0.00949694; 0.0654478]
"pressure EigenMode7" => [-0.995189; -0.991425; … ; 0.954782; 0.884704]
"pressure EigenMode4" => [0.997382; 0.99621; … ; 0.062434; 0.0776241]
The keys are strings and the values are arrays. We can access the values through their keys by using []. For example, to take the third eigenmode:
field_data["pressure EigenMode3"]
This will print on the screen some of the components of the 17949×1 Array{Float64,2} array stored at that key, like this:
17949×1 Array{Float64,2}:
-0.996736720932893
-0.9965775469359124
0.9952550721855146
0.9955974894571784
-0.9962988443020222
-0.9970952870712588
0.9962577599200265
0.9963763219040973
0.9886986673179681
0.9635314855864332
0.9216683604448656
0.8639177119766942
0.7912345455731867
0.7051217070607284
0.6073530288526735
0.4986678330741782
0.38180133279153655
0.2582755976433679
0.1302175588954984
4.819432687419732e-5
-0.13013871030616753
-0.2581526932836525
-0.38163997389416504
⋮
0.5357278504017688
-0.4877328178176225
-0.6894571970644491
-0.4960089540714818
0.5328936014456157
0.5231053593217266
0.5466584719087205
0.33176704659759537
0.5193567081082694
-0.07018174606060222
-0.10554585537059677
-0.05141096021986015
0.22818801520373563
0.3401259380517371
0.42541152891813244
0.3373007444508751
0.22731136783371517
0.2575703033377205
0.3640669793098258
0.44838052883949125
0.4187116628025465
0.4469962968130989
0.32913283350300465
Each scalar value is the value of the scalar field at a node in the mesh. What about the nodes coordinates?
The nodes coordinates are accessible through the attribute points:
coordinates = data.points
The REPL will show you that the points attribute contains a matrix. The number of rows is the same as the number of nodes, while the matrix has 3 columns. It is easy to figure out that each row stores the $x$, $y$ and $z$ values of the mesh point. This is how the matrix will look like in the REPL:
17949×3 Array{Float64,2}:
5.0 4.0 0.0
5.0 4.0 3.0
5.0 0.0 0.0
5.0 0.0 3.0
0.0 4.0 0.0
0.0 4.0 3.0
0.0 0.0 0.0
0.0 0.0 3.0
0.0 0.166667 3.0
0.0 0.333333 3.0
0.0 0.5 3.0
0.0 0.666667 3.0
0.0 0.833333 3.0
0.0 1.0 3.0
0.0 1.16667 3.0
0.0 1.33333 3.0
0.0 1.5 3.0
0.0 1.66667 3.0
0.0 1.83333 3.0
0.0 2.0 3.0
0.0 2.16667 3.0
0.0 2.33333 3.0
0.0 2.5 3.0
⋮
1.24967 1.27832 1.34488
3.05 2.65043 1.31958
2.74866 2.97221 1.32749
2.76474 2.66275 1.5849
2.891 1.28258 1.58422
3.06091 1.29711 1.49163
2.90482 1.26175 1.42878
2.79043 1.56873 1.60513
2.89285 1.30299 1.76927
2.63602 2.08969 1.596
2.6083 2.1349 1.43203
2.81514 2.06576 1.4742
3.19625 1.70616 1.5382
3.11367 1.55699 1.45224
3.13185 1.4389 1.60911
2.05778 1.56065 1.6477
2.07417 1.70751 1.60579
2.05608 1.66748 1.44089
2.0519 1.5242 1.4963
2.05755 1.40648 1.59698
2.08012 1.44836 1.73683
2.27767 1.40842 1.43525
2.12315 1.57164 1.77477
Manipulating the Fields
Having access to the field is enough for a lot of tasks, for example computing metrics, or interpolating along a line, or other things that might be of interest. However, we might wish to manipulate the data and save them somewhere, maybe to do some further processing with ParaView or something like that. That is easy to do, but there are a couple of counter intuitive caveats to be aware of.
Let’s assume we are happy to overwrite the original data. For example, we maybe know that there was an error when computing the fifth eigenmode and all values are wrong by a factor of $10$. So, we maybe think that this will work:
data.point_data["pressure EigenMode5"] = 10 .* data.point_data["pressure EigenMode5"]
However, by printing data.point_data["pressure EigenMode5"] on the REPL before and after running the line above (which you can do by simply typing data.point_data["pressure EigenMode5"]), you will see that the line above does not overwrite the pressure EigenMode5 scalar field with the new data. To overwrite the data we have to make a new dictionary:
new_dict = data.point_data
new_dict["pressure EigenMode5"] = 10 .* new_dict["pressure EigenMode5"]
data.point_data = new_dict
By doing so we update the original data object with all the new scalar fields. We can now save it to file with the write function:
meshio.write("example.vtu", data)
You will be able to load example.vtu in ParaView and observe the corrected scalar field.
But maybe you do not want to overwrite the original data. Maybe you need to retain it for other calculations. So, you might think that doing this will prevent overwriting:
# Let's reload the original data
data = meshio.read("case_t0001.vtu")
# Let's try to copy it in a new variable:
data_modif = data
# And overwrite one of the scalar fields
new_dict = data_modif.point_data
new_dict["pressure EigenMode5"] = 10 .* new_dict["pressure EigenMode5"]
data_modif.point_data = new_dict
By printing data_modif.point_data["pressure EigenMode5"] on the REPL you will see that this worked. However, it had a spooky side effect: also the data in data.point_data["pressure EigenMode5"] were overwritten. Turns out that, if data_modif = data is a copy, it is a shallow one at best, and the two objects have references to the same data. Changing one of them changes them both. Interestingly, deepcopy does not help and results in the same behaviour. What works however is to create data_modif by reading the original data from file again:
data = meshio.read("case_t0001.vtu")
data_modif = meshio.read("case_t0001.vtu")
new_dict = data_modif.point_data
new_dict["pressure EigenMode5"] = 10 .* new_dict["pressure EigenMode5"]
data_modif.point_data = new_dict
Conclusion
In this short episode we very briefly introduced meshio from a Julia perspective. We shown how to load VTU files and manipulate them, as well as a couple of caveats that are important when it comes to data manipulation, overwriting or preventing overwrite of data. This information should be enough to understand the validation code used in some of the previous episodes, as well as to inform the user about the basics to write simple scripts to study and manipulate VTU files programmatically with Julia.
Featured Picture
This article featured picture is from Meshio.
License Information
This work is licensed under a Creative Commons Attribution 4.0 International License.