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], 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)
# Kinematic boundary condition for free surface leads to boundary displacement
# 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()
```