Two ways to use MeshFunction to mark domain boundaries

mesh
1

Hello everyone,

I’m working on a free surface problem, in which the top boundary of a rectangle is a surface. As I calculate the surface curvature, I couldn’t understand why the results were zero, even though the boundary had significant curvature. Suppose I have 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)

Now, I can mark the top boundary in two different ways:

#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

# 2
bmf = MeshFunction("size_t", mesh, d - 1, 0)
bmf.set_all(0)
Top().mark(bmf, 1)

ds = Measure('ds', domain = mesh, subdomain_data = bmf) # for method 2

By using the first method, I can easily use the ALE method to update the boundary coordinate, as the marker works well. But If I use the measure ‘ds’ to integrate over the marked boundary, I get zero. If I use the second method to mark the boundaries, I obtain a non-zero value for the boundary curvature. I deform the boundary using the function:

def init_surface(mesh, bmesh, bmf, L, bflag):
    '''
    Move specified boundary in a senoidal manner.
    '''

    for i in range(bmesh.num_vertices()):
        if (bmf[i] == bflag):
            bmesh.coordinates()[i,1] -= 0.05*np.sin(4*(np.pi/L)*bmesh.coordinates()[i,0])

    ALE.move(mesh, bmesh)
    mesh.bounding_box_tree().build(mesh)

The function init_surface works well when I use the first method to mark the boundaries. However, using the second method I get a deformation on the whole domain, which is not what it is desired.

My question is: what are the differences between those two methods? I’m studying the documentation, but I still haven’t found anwsers. Many thanks in advance.

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
3

Thanks a lot for the fast and detailed reply!