Error defining surface and measures for a 3D mesh imported from Gmsh

General
4

Edge integrals are not supported in 3D meshes without using MeshView.
You have corrected the integral measure for the surfaces (ds, not dS)
Consider the following:

from dolfin import *
import meshio


def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={
                           cell_type: cells}, cell_data={"name_to_read": [cell_data]})
    return out_mesh


######################################################################################################################
msh = meshio.read("mesh.msh")

######################################################################################################################

tetra_mesh = create_mesh(msh, "tetra", True)
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", tetra_mesh)
meshio.write("surface.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)

mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(), "mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("surface.xdmf") as infile:
    infile.read(mvc, "name_to_read")
sf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-2)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

dx_custom = Measure("ds", domain=mesh, subdomain_data=sf)    # Surface measure
dv_custom = Measure("dx", domain=mesh, subdomain_data=cf)    # Volume measure

# Should be 2.5
print(f"Volume = {assemble(Constant(1.0)*dv_custom(19))}")
# Should be 5
print(f"frontsurface = {assemble(Constant(1.0)*dx_custom(13))}")
# Should be 0.25
print(f"leftsurface = {assemble(Constant(1.0)*dx_custom(14))}")
# Should be 0.25
print(f"rightsurface = {assemble(Constant(1.0)*dx_custom(15))}")
# Should be 5
print(f"topsurface = {assemble(Constant(1.0)*dx_custom(16))}")
# Should be 5
print(f"bottomsurface = {assemble(Constant(1.0)*dx_custom(17))}")
# Should be 5
print(f"backsurface = {assemble(Constant(1.0)*dx_custom(18))}")


def assemble_edge(tag):
    mv = MeshView.create(mf, tag)
    return assemble(Constant(1.0)*dx(domain=mv))


# Should be 10
print(f"frontbottomline = {assemble_edge(4)}")
# Should be 10
print(f"fronttopline = {assemble_edge(3)}")
# Should be 0.5
print(f"frontleftline = {assemble_edge(1)}")
# Should be 0.5
print(f"frontrightline = {assemble_edge(2)}")
# Should be 0.5
print(f"lefttopline = {assemble_edge(5)}")
# Should be 0.5
print(f"leftbottomline = {assemble_edge(6)}")
# Should be 0.5
print(f"righttopline = {assemble_edge(7)}")
# Should be 0.5
print(f"rightbottomline = {assemble_edge(8)}")
# Should be 10
print(f"backtopline = {assemble_edge(9)}")
# Should be 10
print(f"backbottomline = {assemble_edge(10)}")
# Should be 0.5
print(f"backleftline = {assemble_edge(11)}")
# Should be 0.5
print(f"backleftline = {assemble_edge(12)}")

returning


Volume = 2.5000000000000027
frontsurface = 5.0
leftsurface = 0.25
rightsurface = 0.25
topsurface = 5.0
bottomsurface = 5.0
backsurface = 5.0
frontbottomline = 10.0
fronttopline = 10.0
frontleftline = 0.5
frontrightline = 0.5
lefttopline = 0.5
leftbottomline = 0.5
righttopline = 0.5
rightbottomline = 0.5
backtopline = 10.0
backbottomline = 10.0
backleftline = 0.5
backleftline = 0.5

For further posts, please remove extra imports and other unused parts of the code.

1 Like