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, 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)
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.