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

[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