Integration over subdomain

4

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_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={
                           "name_to_read": [cell_data]})
    return out_mesh


triangle_mesh = create_mesh(mesh_from_file, "triangle")
meshio.write("facet_mesh.xdmf", triangle_mesh)

tetra_mesh = create_mesh(mesh_from_file, "tetra")
meshio.write("mesh.xdmf", tetra_mesh)

# Read mesh
mesh = Mesh()
mvc_v = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
    infile.read(mvc_v, "name_to_read")
markers = cpp.mesh.MeshFunctionSizet(mesh, mvc_v)

mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
sub_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds_c = Measure("ds", domain=mesh, subdomain_data=sub_domains)
print(assemble(1*ds_c(30)))

yielding

3.1179398656388186
1 Like