# Obtain velocity normal to boundary

Hi everyone. I hope you can guide me to further my studies. I’m working on a 2D free surface problem with a sine wave as initial shape for the top boundary. I solve the 2D Navier-Stokes equations using a mixed formulation along with 2 boundary conditions for the free surface (namely a stress condition “`n.T = sigma*kappa`” and a kinematic condition for the velocity at the interface “`u.n = V.n`”.

I then take the solution of the system (i.e. the velocity vector) and project it onthe surface coordinate system `(normal n, tangential t)`. I calculate the surface displacement with the velocity normal to the surface, in this manner:

``````u_normal = inner(n, u)*n
``````

However the type of `u_normal` is a `ComponentTensor`, and I can’t figure out how to access the indices of this object. Is there a more intelligent way to perform such operation? I attach below a MWE.

``````from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

class Top(SubDomain):
def __init__(self):
SubDomain.__init__(self)
def inside(self, x, on_boundary):
return near(x, 1)

# Arbitrary time step to get displacement
dt = 0.01
# Create mesh and boundary mesh
mesh =  UnitSquareMesh(20,20)
boundary_mesh = BoundaryMesh(mesh, "exterior", True)

# Mesh markers    mesh_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)    boundary_markers = MeshFunction("size_t", boundary_mesh, 0)

# Mark mesh
mesh_markers.set_all(0)
Top().mark(mesh_markers, 1)

# mark boundary_mesh
boundary_markers.set_all(0)
Top().mark(boundary_markers, 1)

# Create function space for bulk dynamics
V = VectorFunctionSpace(mesh, "CG", 2)
u = Function(V)

# Assign random values to bulk velocity    u.vector().set_local(np.random.rand(u.vector().local_size()))
# Create function space for interface dynamics
D= VectorFunctionSpace(boundary_mesh, "CG", 1)
displacement = Function(D)
displacement.set_allow_extrapolation(True)

# Get normal
n = FacetNormal(mesh)

# Obtain velocity normal vector
#! NOT WORKING
u_normal = n * inner(n, u)

# Interpolate mesh values into boundary function
u_boundary = Function(D)
u_boundary.set_allow_extrapolation(True)

# Interpolate velocity normal vector into boundary function
u_boundary.interpolate(u_normal)
# displacement = u_boundary*dt    top_disp = DirichletBC(D, dt*u_boundary, boundary_markers, 1)

# apply DirichletBC into displacement function
top_disp.apply(displacement.vector())

# move mesh
ALE.move(boundary_mesh, displacement)
ALE.move(mesh, boundary_mesh)
mesh.bounding_box_tree().build(mesh)

plot(mesh); plt.show()``````

As the `FacetNormal` as implied by its name, is defined on each Facet, and not for each degree of freedom, you have to approximate the normal at the degrees of freedom (for instance at corners).

For your particular case, the problem is that you are trying to interpolate an `ComponentTensor`, that is an expression using ufl syntax. To work around this, I would suggest using an `UserExpression` as follows:

``````n = FacetNormal(mesh)
# Interpolate mesh values into boundary function
u_boundary = Function(D)
u_boundary.set_allow_extrapolation(True)

# Approximate facet normal in a suitable space using projection
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 2)
u_ = TrialFunction(V)
v_ = TestFunction(V)
a = inner(u_,v_)*ds
l = inner(n, v_)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)

A.ident_zeros()
nh = Function(V)

solve(A, nh.vector(), L)

class normal_u(UserExpression):
def __init__(self, **kwargs):
super().__init__(**kwargs)

def eval(self, values, x):
n_eval = nh(x)
u_eval = u(x)
un = (u_eval*n_eval + u_eval*n_eval)
values = un * n_eval
values = un * n_eval

def value_shape(self):
return (2,)

u_boundary.interpolate(normal_u())
``````