Two ways to use MeshFunction to mark domain boundaries

mesh
2

You first suggested method, as far as I can tell, should return an error, as you are mixing the the domain mesh with a sudomain data from another mesh. For me, running the following:

from dolfin import * # FEniCs library
from mshr import *
import numpy as np

class Top(SubDomain): # This will later on become a free surface
    def __init__(self):
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        return near(x[1], H)

L = 2*np.pi  # Length of channel
H = 1       # Height of channel

# Create mesh
channel = Rectangle(Point(0,0), Point(L, H))
mesh = generate_mesh(channel, 50)

#1
bmesh = BoundaryMesh(mesh, "exterior", True)
bmf = MeshFunction("size_t", bmesh, 0)
bmf.set_all(0)
Top().mark(bmf, 1)

ds = Measure('ds', domain = mesh, subdomain_data = bmf) # for method 1
print(assemble(1*ds))

yields:

*** Error:   Unable to extract mesh from form.
*** Reason:  Non-matching meshes for function spaces and/or measures.
*** Where:   This error was encountered inside Form.cpp.
*** Process: 0

as is expected.
The second method, using:

# 2
bmf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
bmf.set_all(0)
Top().mark(bmf, 1)
ds = Measure('ds', domain = mesh, subdomain_data = bmf) 
print(assemble(1*ds(1)))

yields 2\pi as expected. As you would only like to update the boundary degrees of freedom, you could consider the following:

V= VectorFunctionSpace(mesh, "CG", 1)
u = Function(V)
bc = DirichletBC(V,  Expression(("0", "0.05*sin(4*(pi/L)*x[0])"), L=L, degree=1), bmf, 1)
bc.apply(u.vector())
ALE.move(mesh, u)
File("mesh.pvd") << mesh
1 Like