Working with Imported Meshes

[stefano.tronci]

Hi! My problem is perhaps very simple, but I could not find an answer yet.

Basically, I am making complex meshes with Salome. These are meshes with various solids inside. To make things simple, I made a simple mesh made of 2 concentric spheres to start things off. Basically, I aim to run acoustics simulations buy having the outer shell acting as an adiabatic absorber.

I pretty much followed Dokken’s excellent tutorial and I came up with this code:

import pathlib
import meshio
import tempfile
import fenics


# This is a 3D mesh with concentring spheres, differently marked.
mesh_path = pathlib.Path('ConcentricSpheres.med')

# We read te MED file with meshio
raw_mesh_data = meshio.read(mesh_path)

# This list will contain FEniCS meshes, in this order: line (1D), triangle (2D) and tetra (3D).
# To make the FEniCS meshes we crate meshio meshes for each entity, save it to a temporary location as xdmf, and then
# we read it back.
meshes = []
with tempfile.TemporaryDirectory() as xdmf_dir:
    target_xdmf_dir = pathlib.Path(xdmf_dir)

    for entity in ['line', 'triangle', 'tetra']:
        cells = raw_mesh_data.get_cells_type(entity)
        cell_data = raw_mesh_data.get_cell_data('cell_tags', entity)
        points = raw_mesh_data.points  # 3D Mesh: No pruning needed.
        edited_mesh_data = meshio.Mesh(points=points, cells={entity: cells}, cell_data={'marker': [cell_data]})

        tmp_path = target_xdmf_dir.joinpath(entity + '_mesh.xdmf')
        meshio.write(tmp_path.absolute().as_posix(), edited_mesh_data)
        mesh = fenics.Mesh()
        with fenics.XDMFFile(tmp_path.absolute().as_posix()) as xdmf_file:
            xdmf_file.read(mesh)
        meshes.append(mesh)

# If everything worked, this will print something like:
# [<dolfin.cpp.mesh.Mesh object at 0x7f9680b05ca8>, <dolfin.cpp.mesh.Mesh object at 0x7f9680b058e0>, <dolfin.cpp.mesh.Mesh object at 0x7f9680b05a98>]
print(meshes)

The mesh to run this code is available here. Things look alright: I will get a dolfin.cpp.mesh.Mesh object for each type of entity (line, triangle and tetra). I can take a look to the xdmf files with ParaView and they seem OK:

However, I am not sure how to make use of my 3 dolfin.cpp.mesh.Mesh objects. So far, I only used a single mesh object, for example:

mesh = fenics.BoxMesh(
        fenics.Point(0, 0, 0),
        fenics.Point(1, 1, 1),
        10,
        10,
        10
    )

Should I merge all the 3 meshes together in one single mesh? If so, how? Or should I just use the tetra one?

P.S.
The various ID numbers coming from the MED file are all negative. I might have read in some other forum post (which I cannot dig out) that this can cause issues. Should I make them all positive?

Cheers!

[dokken]

Instead of reading in the 1D and 2D meshes as meshes, they should be considered as MeshValueCollections. See for instance

Integration over subdomain

The boundary with tag 30 in your mesh is not an interior surface. It is an external surface, as there are only cells on one side of the facets. Thus you can use the ds Measure, without the “+” and “-” restriction to integrate over the surface. For instance consider the minimal code (using your mesh): from dolfin import * import meshio # Convert mesh mesh_from_file = meshio.read("sim_2.msh") def create_mesh(mesh, cell_type, prune_z=False): cells = mesh.get_cells_type(cell_type) cell…

[stefano.tronci]

Thank you very much for your help @dokken !

I followed your example, with a few changes:

• I used int instead of size_t as all of my IDs are signed integers;
• I used fenics.MeshFunction instead of dolfin.cpp.mesh.MeshFunctionInt, mainly for consistency.

Do you have any comment about this?

The code runs. I assume that now I will be able to setup any model by using the mesh, the markers and the subdomains as I see fit?

My update code is below, for reference.

import pathlib
import meshio
import tempfile
import fenics
# import dolfin.cpp.mesh


# This is a 3D mesh with concentring spheres, differently marked.
mesh_path = pathlib.Path('ConcentricSpheres.med')

# We read te MED file with meshio
raw_mesh_data = meshio.read(mesh_path)

# This list will contain FEniCS meshes, in this order: line (1D), triangle (2D) and tetra (3D).
# To make the FEniCS meshes we crate meshio meshes for each entity, save it to a temporary location as xdmf, and then
# we read it back.
mesh = fenics.Mesh()
mvc_3d = fenics.MeshValueCollection('int', mesh, 3)
mvc_2d = fenics.MeshValueCollection('int', mesh, 2)
mvc_1d = fenics.MeshValueCollection('int', mesh, 1)
with tempfile.TemporaryDirectory() as xdmf_dir:
    target_xdmf_dir = pathlib.Path(xdmf_dir)

    for entity in ['tetra', 'triangle', 'line']:
        cells = raw_mesh_data.get_cells_type(entity)
        cell_data = raw_mesh_data.get_cell_data('cell_tags', entity)
        points = raw_mesh_data.points  # 3D Mesh: No pruning needed.
        edited_mesh_data = meshio.Mesh(points=points, cells={entity: cells}, cell_data={'marker': [cell_data]})

        tmp_path = target_xdmf_dir.joinpath(entity + '_mesh.xdmf')
        meshio.write(tmp_path.absolute().as_posix(), edited_mesh_data)
        with fenics.XDMFFile(tmp_path.absolute().as_posix()) as xdmf_file:
            if entity == 'tetra':
                xdmf_file.read(mesh)
                xdmf_file.read(mvc_3d, 'marker')
            if entity == 'triangle':
                xdmf_file.read(mvc_2d, 'marker')
            if entity == 'line':
                xdmf_file.read(mvc_1d, 'marker')

# markers = dolfin.cpp.mesh.MeshFunctionInt(mesh, mvc_3d)
markers = fenics.MeshFunction('int', mesh, mvc_3d)
# sub_domains_2d = dolfin.cpp.mesh.MeshFunctionInt(mesh, mvc_2d)
sub_domains_2d = fenics.MeshFunction('int', mesh, mvc_2d)
# sub_domains_1d = dolfin.cpp.mesh.MeshFunctionInt(mesh, mvc_1d)
sub_domains_1d = fenics.MeshFunction('int', mesh, mvc_1d)

ds_c = fenics.Measure('ds', domain=mesh, subdomain_data=sub_domains_2d)
print(ds_c(-9))

[dokken]

@stefano.tronci

Do you have any comment about this?

This looks fine to me. Good luck with your simulations.

[stefano.tronci]

@stefano.tronci

P.S.
The various ID numbers coming from the MED file are all negative. I might have read in some other forum post (which I cannot dig out) that this can cause issues. Should I make them all positive?

Hi! I found that negative ID numbers are indeed an issue, as it will make the compiler make names with - in it, that then cause error. It is easy to change all IDs to be positive though. I just applied the same offset to them all, and used size_t as a type. Updated code below.

import pathlib
import meshio
import tempfile
import numpy as np
import fenics
# import dolfin.cpp.mesh


# This is a 3D mesh with concentring spheres, differently marked.
mesh_path = pathlib.Path('ConcentricSpheres.med')

# We read te MED file with meshio
raw_mesh_data = meshio.read(mesh_path)
cell_tag_offset = np.abs(np.min(list(raw_mesh_data.cell_tags.keys())))

# This list will contain FEniCS meshes, in this order: line (1D), triangle (2D) and tetra (3D).
# To make the FEniCS meshes we crate meshio meshes for each entity, save it to a temporary location as xdmf, and then
# we read it back.
mesh = fenics.Mesh()
mvc_3d = fenics.MeshValueCollection('size_t', mesh, 3)
mvc_2d = fenics.MeshValueCollection('size_t', mesh, 2)
mvc_1d = fenics.MeshValueCollection('size_t', mesh, 1)
with tempfile.TemporaryDirectory() as xdmf_dir:
    target_xdmf_dir = pathlib.Path(xdmf_dir)

    for entity in ['tetra', 'triangle', 'line']:
        cells = raw_mesh_data.get_cells_type(entity)
        cell_data = raw_mesh_data.get_cell_data('cell_tags', entity)
        points = raw_mesh_data.points  # 3D Mesh: No pruning needed.
        edited_mesh_data = meshio.Mesh(
            points=points,
            cells={entity: cells},
            cell_data={'marker': [cell_data + cell_tag_offset]}
        )

        tmp_path = target_xdmf_dir.joinpath(entity + '_mesh.xdmf')
        meshio.write(tmp_path.absolute().as_posix(), edited_mesh_data)
        with fenics.XDMFFile(tmp_path.absolute().as_posix()) as xdmf_file:
            if entity == 'tetra':
                xdmf_file.read(mesh)
                xdmf_file.read(mvc_3d, 'marker')
            if entity == 'triangle':
                xdmf_file.read(mvc_2d, 'marker')
            if entity == 'line':
                xdmf_file.read(mvc_1d, 'marker')

markers = fenics.MeshFunction('size_t', mesh, mvc_3d)
sub_domains_2d = fenics.MeshFunction('size_t', mesh, mvc_2d)
sub_domains_1d = fenics.MeshFunction('size_t', mesh, mvc_1d)

ds_c = fenics.Measure('ds', domain=mesh, subdomain_data=sub_domains_2d)
print(ds_c(9))

After that, it seems I can at least solve my prototype sims (now I will check for accuracy).