Wrong result of inner product evaluated at boundaries

variational formulation
1

Dear all,

I am having a problem when computing an inner product over the domain boundaries.
My code runs without errors, but I think that the output is wrong because I am combining functions that belong to different finite element spaces.

The definition of the function space and solution vector is:

<\boldsymbol{\sigma}, \boldsymbol{\mathcal{U}}>_{\Gamma} = \int_{\Gamma} \left[\left(p\hat{\boldsymbol{n}} + \frac{1}{Re}\hat{\boldsymbol{n}} \cdot \nabla \boldsymbol{u} \right) \cdot \boldsymbol{\mathcal{U}}\right] d\Gamma
where p is pressure, \boldsymbol{u} is velocity, \hat{\boldsymbol{n}} is the unit normal vector to the boundaries, \boldsymbol{\mathcal{U}} is the velocity at the boundaries and \Gamma defines the boundaries.

# Function Space
el1 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
el2 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el1, el2])
Y = fem.FunctionSpace(domain, mel)
# Solution vector
q_n = fem.Function(Y)
u_n = q_n.sub(0)
p_n = q_n.sub(1)

Then, I compute the inner product as follows:

Fun = fem.form( ufl.inner((p_n*n + (1/Re)*ufl.dot(n,ufl.grad(u_n))), v)*ufl.ds )
Fun_assemble = MPI.COMM_WORLD.gather(fem.assemble_scalar(Fun), root=0)

As I said, my belief is that I’m not projecting p_n into the function space of u_n, but it might also be that I’m not expressing the inner product properly.
Does anybody know what I am doing wrong?

Please, find below a minimal example:

from dolfinx import fem, mesh, io, plot, cpp, graph
import ufl
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
rank = MPI.COMM_WORLD.rank

# Define mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0, 0],[1, 0.1]], [100, 10])
n = ufl.FacetNormal(domain)

# Define function space
el1 = ufl.VectorElement("Lagrange", domain.ufl_cell(), 2)
el2 = ufl.FiniteElement("Lagrange", domain.ufl_cell(), 1)
mel = ufl.MixedElement([el1, el2])
Y = fem.FunctionSpace(domain, mel)

# Define boundaries
def boundary_function_L(x):
    return np.isclose(x[0], 0)
def boundary_function_R(x):
    return np.isclose(x[0], 1)
def boundary_function_B(x):
    return np.isclose(x[1], 0)
def boundary_function_T(x):
    return np.isclose(x[1], 0.1)

# Define solution vector
q_n = fem.Function(Y)
u_n = q_n.sub(0)
p_n = q_n.sub(1)

# Define trial and test functions
u, p = ufl.TrialFunctions(Y)
v, q = ufl.TestFunctions(Y)

# Inner product
Re = 100
Fun = fem.form( ufl.inner((p_n*n + (1/Re)*ufl.dot(n,ufl.grad(u_n))), v)*ufl.ds )
Fun_assemble = MPI.COMM_WORLD.gather(fem.assemble_scalar(Fun), root=0)

Note that here Fun_assemble is going to be zero because u_n and p_n have been just initialised, but in my full code it’s also zero even though u_n and p_n are different from zero.

Many thanks in advance.

2

This is not a scalar form, as there is a test_function in the form. You should use fem.petsc.assemble_vector

3

Thanks for the quick response, @dokken !
But, isn’t the following expression a scalar?

<\boldsymbol{\sigma},\boldsymbol{\mathcal{U}}>_{\Gamma} = \int_{\Gamma}\left[\left(p\hat{\boldsymbol{n}}+\frac{1}{Re}\hat{\boldsymbol{n}}\cdot\nabla\boldsymbol{u}\right)\cdot\boldsymbol{\mathcal{U}}\right]d\Gamma

4

For each input \mathcal{U} this expression will generate a scalar. However, each U is a function in a given function space, which you would take the inner product with the result of fem.petsc.assemble_vector