Vertically integrate the norm of a 3D Vector Function.

Question

I have a vector-valued Function $\vec{u}$ defined on a 3D grid with subdomains marked on the boundary (3 on the base). I would like to vertically integrate the norm of these vectors from the bed to the surface. For this, I solve:

$$ v = \int_b^s ||\vec{u}|| dz, $$
$$ \implies \frac{\partial v}{\partial z} = ||\vec{u}||, \hspace{5mm} v(b) = 0, $$

where $b$ is the height of the base and $s$ is the height of the surface.

However, I am getting the value of the function $u$ on the surface, not the integral which I require.

from dolfin import *

mesh = UnitCubeMesh(10,10,10)
ff   = FacetFunction("size_t", mesh, 0)

# calculate boundaries :
for f in facets(mesh):
  n   = f.normal()
  tol = 1e-3
  # surface :
  if n.z() >= tol and f.exterior():
    ff[f] = 1
  # base :
  elif n.z() <= -tol and f.exterior():
    ff[f] = 2

# define function space :
Q   = FunctionSpace(mesh, "CG", 1)

# a constant throughout the domain for testing :
one = Constant(1)

def vert_integrate(u, ff, Q):
  phi    = TestFunction(Q)
  v      = TrialFunction(Q)
  bc     = DirichletBC(Q, 0.0, ff, 2)
  a      = v.dx(2) * phi * dx
  L      = u * phi * dx
  v      = Function(Q)
  solve(a == L, v, bc)
  return v

# plot solution :
plot(vert_integrate(one, ff, Q))
interactive()

Can anyone help?

Comments

Can you post a short and complete code demonstrating the issue?

---

updated above, thanks.

Answer 1

Your answer is actually correct, the integral is actually equal to one:

$$ \int_0^1 1 dz = [z]_0^1 = 1 - 0 = 1. $$

Note that if you make the domain larger in the $z$-direction, you will have the proper integral as well :

$$ \int_0^{10} 2 dz = [2z]_0^{10} = 2 \cdot 10 - 0 = 20, $$

as demonstrated in the code snippet below:

from dolfin import *

mesh = BoxMesh(0, 0, 0, 1, 1, 20, 10, 10, 10)
ff   = FacetFunction("size_t", mesh, 0)

# calculate boundaries :
for f in facets(mesh):
  n   = f.normal()
  tol = 1e-3 
  # surface :
  if n.z() >= tol and f.exterior():
    ff[f] = 1
  # base :
  elif n.z() <= -tol and f.exterior():
    ff[f] = 2

# define function space :
Q   = FunctionSpace(mesh, "CG", 1)

def vert_integrate(u, ff, Q):
  phi    = TestFunction(Q)
  v      = TrialFunction(Q)
  bc     = DirichletBC(Q, 0, ff, 2)
  a      = v.dx(2) * phi * dx
  L      = u * phi * dx
  v      = Function(Q)
  solve(a == L, v, bc)
  return v

# plot solution :
plot(vert_integrate(Constant(2), ff, Q))
interactive()

Comments

whoa, thanks! How clumsy of me!