Loading meshes

[bhaveshshrimali]

Hi @adtzlr ,

Is it possible to load an external mesh via meshio (I do see that here) but not sure, e.g.

from meshio import xdmf
import felupe as fem

meshio_mesh = xdmf.read("trialMesh.xdmf")
mesh = fem.Mesh(
    points=meshio_mesh.points,
    cells=meshio_mesh.get_cells_type("tetra"),
    cell_type=meshio_mesh.cells[0].type,
)

However this complains the cell volumes being less than zero for some of the cells, I tried flipping as suggested, namely,

region_u = fem.RegionQuadraticTetra(mesh=mesh)
mesh = mesh.flip(np.any(region_u.dV < 0, axis=0))
region_u = fem.RegionQuadraticTetra(mesh=mesh)
region_p = fem.RegionTetra(
    mesh=fem.mesh.dual(mesh, disconnect=False),
    quadrature=region_u.quadrature,
    grad=False,
)

but this seems to just invert the ones having positive volume to negative. Would you know how to properly load an external mesh, or could point me to an example

For a test mesh I used:

from skfem import MeshTet
mesh = MeshTet().refined(2)
mesh.save("trialMesh.xdmf")

[adtzlr]

Hi, I created an issue #624. This is fixed in #625.

from skfem import MeshTet

mesh_skfem = MeshTet().refined(2)
mesh_skfem.save("trialMesh.xdmf")

from meshio import xdmf
import felupe as fem
import numpy as np

meshio_mesh = xdmf.read("trialMesh.xdmf")
mesh = fem.Mesh(
    points=meshio_mesh.points,
    cells=meshio_mesh.get_cells_type("tetra"),
    cell_type=meshio_mesh.cells[0].type,
)

# ensure positive cell volumes (before adding midpoints on edges!)
region = fem.RegionTetra(mesh)
mesh = mesh.flip(np.any(region.dV < 0, axis=0))
region.reload(mesh)
assert np.all(region.dV > 0)

region_u = fem.RegionQuadraticTetra(mesh.add_midpoints_edges())
region_p = fem.RegionTetra(
    mesh=fem.mesh.dual(mesh, disconnect=False),
    quadrature=region_u.quadrature,
    grad=False,
)

P.S.: It is also possible to use container = mesh.read(filename) and select the desired mesh from the list of meshes of the mesh container.

>>> import felupe as fem

>>> container = fem.mesh.read(filename="trialMesh.xdmf")
>>> container.meshes
[<felupe Mesh object>
   Number of points: 115
   Number of cells:
     tetra: 320]

>>> container.meshes[0].imshow()

[bhaveshshrimali]

wonderful, thanks a lot

[adtzlr]

I just pushed v7.18.0 to PyPI with the fix included.