# Function in boundary condition, defined on boundary only

Greetings. I have some function. It is included in the right part of the system of nonlinear PDEs under the normal derivative. The problem is that I have a function defined on boundary only. And multiplication on the test function gives me `TypeError: '<' not supported between instances of 'Mesh' and 'Mesh'`

``````import itertools

from utilities import Normal

# noinspection PyUnresolvedReferences
from dolfin import dx, ds
from dolfin import *

class Normal(UserExpression):
def __init__(self, mesh, dimension, **kwargs):
self.mesh = mesh
self.dimension = dimension
super().__init__(**kwargs)

def eval_cell(self, values, x, ufc_cell):
values[0] = 0
values[1] = 0
for i in range(self.dimension):

if x[i] < DOLFIN_EPS:
values[i] = -1
if x[i] + DOLFIN_EPS > 1:
values[i] = 1
for i, j in itertools.combinations(list(range(self.dimension)), 2):
if abs(values[i]) == abs(values[j]) == 1:
values[0] *= 0.5
values[1] *= 0.5

def value_shape(self):
return (2,) if self.dimension == 2 else (3,)

omega = UnitSquareMesh(10, 10)
omega_b = BoundaryMesh(omega, 'exterior')
finite_element = FiniteElement("Lagrange", omega.ufl_cell(), 2)

state_space = FunctionSpace(omega, finite_element * finite_element)
simple_space = FunctionSpace(omega, finite_element)
vector_space = VectorFunctionSpace(omega, 'Lagrange', 2)
boundary_vector_space = VectorFunctionSpace(omega_b, 'Lagrange', 1)
boundary_simple_space = FunctionSpace(omega_b, 'Lagrange', 1)
normal = interpolate(Normal(omega, dimension=2), boundary_vector_space)
normal_2 = FacetNormal(omega_b)
theta = project(
Expression('pow((x[0]-0.5), 2) - 0.5*x[1] + 0.75', element=simple_space.ufl_element()),
simple_space
)

v, h = TestFunctions(state_space)
state = Function(state_space)
a, b = split(state)

a_eq = inner(grad(a), grad(v)) * dx + inner(a ** 4 - b, v) * dx
b_eq = inner(grad(b), grad(h)) * dx + inner(b - a ** 4, h) * dx
a_src = 0.3 * v * ds
b_src = theta_n * h * ds
solve(a_eq + b_eq - a_src - b_src == 0, state)
``````

What am I missing? How can I include boundary-only function into the boundary conditions?

You cannot mix meshes in a single assembly, one of the reasons is that the integration measure is not well defined in that case.
Could you please explain why you cannot use

``````n = FacetNormal(omega)
b_src = theta_n * h * ds
``````
1 Like

To your question - I have a few reasons:

1. Issue not with the normal derivative itself but rather with â€śfunction defined on boundary onlyâ€ť - on my actual case itâ€™s a 3-rd type boundary condition. So I have to use a combination of the normal derivative plus some other function. Yeah, I guess I could refactor my code a bit and maybe get a result. I will do and report what happened.

2. I would like to be able to manipulate freely with my functions to trace whatâ€™s going on. `FacetNormal` gives `Integral of type cell cannot contain a ReferenceNormal` outside of the problem formulation scope (Iâ€™m sorry if Iâ€™m wrong here) so I personally want to avoid to use it.

I have previously posted a recipe on how to project the FacetNormal to an appropriate volume function space, see: How to plot normal unit vector of faces in a 2D mesh? - #2 by dokken
This function would only give you a normal on the boundary, 0 at all other dofs.

I have to check this out. Thank you, sir.

It works amazingly nice. Finally I could operate with normal derivative in the same way as any other function.

``````def normal_for_square_mesh():
mesh = UnitSquareMesh(100, 100)
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
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)
File("nh.xml") << nh
return nh

normal_function = []

def get_normal_derivative(function):
if not normal_function:
normal_function.append(normal_for_square_mesh())
normal = normal_function[0]
mesh = UnitSquareMesh(100, 100)
f = FunctionSpace(mesh, "CG", 1)
function = interpolate(function, f)
V = VectorFunctionSpace(mesh, "CG", 1)
return project(inner(normal, project(grad(project(function, f)), V)), f)
``````