From .vtk to .xdmf. Reading mesh infos and setting up boundaries

[Evandro_Gaio]

Hi everyone. I’m trying to validate my computational model using a .vtk mesh.
This .vtk file has everything I need, i.e:

So I used mesh.io to convert this .vtk file. I tried for .xdmf and .XML formats. I’ve always used.XML format, but in the conversion, it says that this format is legacy, and should I consider using .xdmf instead.

The 1st conversion, .vtk to . XML, I got 4 files: mesh, fibers, ID, and sheets, as I expected. But when I tried to read the files like this:

mesh = Mesh(“mesh-1.xml”)
ids = MeshFunction(“size_t”, mesh, “mesh-1_ID.xml”)
fibers = MeshFunction(“size_t”, mesh, “mesh-fibers.xml”)
sheets = MeshFunction(“size_t”, mesh, “mesh-sheets.xml”)

besides the mesh, I got the following error:

Error: Unable to read mesh function from XML file
Reason: Type mismatch reading XML MeshFunction. MeshFunction type is “unit”, but file type is “int” (for ids) and “float” (for fibers/sheets). I even tried to change “size_t” to int or float respectively but didn’t work.

So I tried the .xdmf format. I got only one file out from the conversion: mesh.vtk → mesh.xdmf. So my question is: How to read this only mesh.xdmf file and get all properties showed in the figure (ids, fibers and sheets) ?

The last question is: I have this scalar field “rho”, the goes from 0 to 1 in the rho direction. I’d like to know how to setup my boundaries for rho=0 and rho=1. How can I use this scalar field to define boundaries?

Any help would be appreciated. Thanks a lot!

[dokken]

I would suggest using meshio, as shown in

Accessing and marking imported boundaries

As you are using a stp file I do not have access to, I cannot really do much more than showing you a minimal working example. Here is a minimal geo file (taken from the gmsh tutorial, called tutorial_t1.geo lc = 1e-1; Point(1) = {0, 0, 0, lc}; Point(2) = {.1, 0, 0, lc}; Point(3) = {.1, .3, 0, lc}; Point(4) = {0, .3, 0, lc}; Line(1) = {1, 2}; Line(2) = {3, 2}; Line(3) = {3, 4}; Line(4) = {4, 1}; Curve Loop(1) = {4, 1, -2, 3}; Plane Surface(1) = {1}; Physical Curve(666) = {1, 2, 4}; Physical S…

To be slightly more specific, read in your vtk mesh (and data ) with meshio.read("filename.vtk"), then create a mesh with each of the data

    mesh = meshio.read("filename.vtk")
    cells = mesh.get_cells_type(cell_type)
    # Code for extracting fibers from the read in mesh goes here 
    strain_data = ....
    fibers_data = ....
    sheets_data = ...
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"strain": [strain_data], "fibers":[fibers_data], ....]})
    if prune_z:
        out_mesh.prune_z_0()
   out_mesh.write("mesh.xdmf")  

which can be read into dolfin as

mesh = Mesh()

with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
    infile.read(mvc, "strain")
strain_func = cpp.mesh.MeshFunctionSizet(mesh, mvc)  # boundary_parts

The code above is pseudo-code, as I do not use the vtk format, and do not have any minimal examples at hand to test the code. Therefore, this might require modification. Note that there are many posts about how to read in meshes and external data on this forum.

For your last question, how is your rho field defined? is it defined as a function of x,y,z? if so, use subdomain on the following form:

 class rho0_boundary(SubDomain):
        def inside(self, x, on_boundary):
            return on_boundary and near(rho(x), 0)

and use the subdomain as shown in: https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=measure

[Evandro_Gaio]

First of all thanks for replaying Dokken. I opened my .xdmf file and it still having the same data from .vtk file, so this is good so far. Then I skipped to the second part, to read the .xdmf file, that is:

<Xdmf Version="3.0"><Domain><Grid Name="Grid"><Geometry GeometryType="XYZ"><DataItem DataType="Float" Dimensions="171860 3" Format="HDF" Precision="4">mesh-1.h5:/data0</DataItem></Geometry><Topology TopologyType="Tetrahedron" NumberOfElements="887489" NodesPerElement="4"><DataItem DataType="Int" Dimensions="887489 4" Format="HDF" Precision="4">mesh-1.h5:/data1</DataItem></Topology><Attribute Name="/data/KCL/Eidolon/P150/P150-0.8mm-vol-final.uvc/UVC/COORDS_PHI.dat" AttributeType="Scalar" Center="Node"><DataItem DataType="Float" Dimensions="171860 1" Format="HDF" Precision="4">mesh-1.h5:/data2</DataItem></Attribute><Attribute Name="/data/KCL/Eidolon/P150/P150-0.8mm-vol-final.uvc/UVC/COORDS_RHO.dat" AttributeType="Scalar" Center="Node"><DataItem DataType="Float" Dimensions="171860 1" Format="HDF" Precision="4">mesh-1.h5:/data3</DataItem></Attribute><Attribute Name="/data/KCL/Eidolon/P150/P150-0.8mm-vol-final.uvc/UVC/COORDS_Z.dat" AttributeType="Scalar" Center="Node"><DataItem DataType="Float" Dimensions="171860 1" Format="HDF" Precision="4">mesh-1.h5:/data4</DataItem></Attribute><Attribute Name="ID" AttributeType="Scalar" Center="Cell"><DataItem DataType="Int" Dimensions="887489 1" Format="HDF" Precision="4">mesh-1.h5:/data5</DataItem></Attribute><Attribute Name="fibres" AttributeType="Vector" Center="Cell"><DataItem DataType="Float" Dimensions="887489 3" Format="HDF" Precision="4">mesh-1.h5:/data6</DataItem></Attribute><Attribute Name="sheets" AttributeType="Vector" Center="Cell"><DataItem DataType="Float" Dimensions="887489 3" Format="HDF" Precision="4">mesh-1.h5:/data7</DataItem></Attribute></Grid></Domain></Xdmf>

so I used:

mesh = Mesh()
 
 with XDMFFile("mesh.xdmf") as infile:
     infile.read(mesh)
     mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
     infile.read(mvc, "/data/KCL/Eidolon/P150/P150-0.8mm-vol-final.uvc/UVC/COORDS_PHI")
strain_func = cpp.mesh.MeshFunctionSizet(mesh, mvc)  # boundary_parts

and worked for “PHI”. However I tried for the others data, like RHO, Z, ID, Fibers and Sheets and got the following error:

*** Error: Unable to recognise cell type.
*** Reason: Unknown value “”.
*** Where: This error was encountered inside XDMFFile.cpp.

Do you have any idea how could it worked for PHI but it doesn’t for the others?

About the RHO to use as a boundary, RHO is a scalar: 0<= RHO <=1.

[dokken]

It is quite hard to help you as I do not have any VTK file with similar data. Could you make a small VTK file with like 4 cells and similar data and post it here (with 3x` encapsulation)? If you do that I’ll have a go at it.

So is rho defined for each facet or each cell of the mesh?

[Evandro_Gaio]

So the mesh file I was trying to use is open access. from:

https://librarysearch.kcl.ac.uk/primo-explore/fulldisplay?vid=kings_v2&docid=44KCL_ALEPH_RDM_DS061309465&context=L

The mesh I trying to use is the 1st one: link

https://emckclac.sharepoint.com/sites/rdm01/rdmdataset/Forms/AllItems.aspx?id=%2Fsites%2Frdm01%2Frdmdataset%2FA%20Virtual%20Cohort%20of%20Twenty-four%20Left-ventricular%20Models%20of%20Ischemic%20Cardiomyopathy%20Patients%2F1-mesh-uvc.vtk&parent=%2Fsites%2Frdm01%2Frdmdataset%2FA%20Virtual%20Cohort%20of%20Twenty-four%20Left-ventricular%20Models%20of%20Ischemic%20Cardiomyopathy%20Patients&p=true&originalPath=aHR0cHM6Ly9lbWNrY2xhYy5zaGFyZXBvaW50LmNvbS86dTovcy9yZG0wMS9FUW9nUkpMamlQMUNsbGJiNUx6V1hQQUJwR0kyQjJZZ0FyTWRMNGs2UDd3NnpnP3J0aW1lPVhySVNOcTgxMlVn

Rho is defined on each point.

[dokken]

The second link does not work, however, I got the mesh from the first link

However, looking at the data, The mesh format isn’t very suitable for xdmf. What I would consider doing is to convert the data to xml, read it into dolfin, then save it as xdmf using dolfin functionality. Especially the cell data that is three dimensional is quite clunky to convert with meshio.

[Evandro_Gaio]

Got it. Thank you for your support anyway. I’ll keep trying to make it works.