Integration over subdomain

1

Hello every body, is there any way to compute surface integrals over interior surface? For example if I had a sphere inside a cube, how it is possible to consider only the sphere surface instead of dx which consider the entire volume?

2

You should use dS with an appropriate mesh function and a volume marker to have a consistent restriction of the “+” and “-”. See: Integrating over an interior surface - #4 by MiroK

3

Maybe is better if Iupload some code, becouse I get 0 as integral when I follow your example.

# Parameters
mu = 1.49
rho = 1259
nu = mu/rho

# 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)

# Build function spaces on Mini element
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
B = FiniteElement("Bubble",   mesh.ufl_cell(), mesh.topology().dim() + 1)
V = VectorElement(NodalEnrichedElement(P1, B))
Q = P1
W = FunctionSpace(mesh, V * Q)

# Define variational problem
(v, q) = TestFunctions(W)
(u, p) = TrialFunctions(W)

# Define boundary conditions
v_inlet = Constant((0, 0, 1))
v_outlet = Constant((0, 0, 0))
p_inlet = Constant(5)
p_outlet = Constant(4)
noslip = Constant((0, 0, 0))

# Define velocity profile
rotation_p1 = ('.5*(x[1]-0.5)', '.5*(x[0]-0)', '0')
rotation_p2 = ('-.5*(x[1]+0.5)', '-.5*(x[0]-0)', '0')
rotation_p3 = ('.5*(x[1]-0)', '.5*(x[0]+0.5)', '0')
rotation_p4 = ('-.5*(x[1]+0)', '-.5*(x[0]-0.5)', '0')
rotation_hull = ('0.7*(x[1]-0.5)', '0.7*(x[1]-0.5)', '0.7*(x[1]-0.5)')

# Define boundary conditions
#bcu_hull = DirichletBC(W.sub(0), Expression(rotation_hull, degree=1), sub_domains, 70)
bcu_inflow = DirichletBC(W.sub(0), v_inlet, sub_domains, 27)
bcp_outflow = DirichletBC(W.sub(1), p_outlet, sub_domains, 29)
bcu_walls = DirichletBC(W.sub(0), noslip, sub_domains, 28)
#bcu_hull = DirichletBC(W.sub(0), Expression(rotation_hull, degree=1), sub_domains, 30)
bcu_hull = DirichletBC(W.sub(0), noslip, sub_domains, 30)

bcu_propeller1 = DirichletBC(W.sub(0), Expression(rotation_p1, degree=1), sub_domains, 31)
bcu_propeller4 = DirichletBC(W.sub(0), Expression(rotation_p4, degree=1), sub_domains, 34)
bcu_propeller2 = DirichletBC(W.sub(0), noslip, sub_domains, 32)
bcu_propeller3 = DirichletBC(W.sub(0), noslip, sub_domains, 33)

bcs = [bcu_inflow, bcu_hull, bcu_propeller1,bcu_propeller2,bcu_propeller3 ,bcu_propeller4, bcu_walls, bcp_outflow]

# Define linear and bilinear forms with stabilized terms
f = Constant((0, 0, 0))
h = 1e-3
beta  = 0.2
delta = beta*h*h
a = (inner(grad(v), grad(u)) - div(v)*p + q*div(u) + \
    delta*inner(grad(q), grad(p)))*dx
L = inner(v + delta*grad(q), f)*dx

# Compute solution
w = Function(W)
solve(a == L, w, bcs, solver_parameters= {'linear_solver': 'mumps'})
#solve(a == L, w, bcs)
u, p = w.split()

def symgrad(u):
    return sym(grad(u))

dS = Measure('dS', domain=mesh, subdomain_data=sub_domains)
K = np.zeros((3,3))
for i in range(3):
    for j in range(3):
        K[i,j] = 2*mu*assemble(inner(symgrad(u('+'))[:,i], symgrad(u('+'))[:,j])*dS(30)) + \
                 2*mu*assemble(inner(symgrad(u('-'))[:, i], symgrad(u('-'))[:, j]) * dS(30))
        print("K_ij =", K[i,j])

What i want to compute is the integral of the sym grad of the solution over the boundaries of the hull which is identified by the identification number: 30. The mesh is findable at this public folder:
https://drive.google.com/drive/folders/1_pOWo9-23KmrvvI6uviZc6Hwmsn-pCQT?usp=sharing

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