### How to solve Poisson equation with tangential boundary condition?

I encounter a problem, which can be expressed as follows.

I want to solve a PDE, briefly, expressed by a Poisson equation.

$\Delta u=f$Δ

`u`=

`ƒ`with boundary condition $u=\left(\nabla\phi\cdot\tau\right)\tau$

`u`=(∇

`ϕ`·

`τ`)

`τ`

Where $\tau$

`τ`is the tangential direction of boundaries, and $\phi$

`ϕ`is a known function or can be solved by the previous step.

I have browse the FEniCS forum for inspire, some similar topics listed bellow,

https://fenicsproject.org/qa/9926/boundary-conditions-in-normal-and-tangent-directions?show=9926#q9926

https://fenicsproject.org/qa/8590/plot-facetnormal?show=8594#a8594

https://fenicsproject.org/qa/6635/possible-error-assembling-lagrange-multipliers-boundaries?show=6635#q6635

https://fenicsproject.org/qa/6177/how-to-compute-point-wise-quantity-along-a-boundary

However, my problem still be remained. Do I miss some important information of the forum topics ?

Could you give me some comment for this problem?

Thank so much!

Best, Hamilton

### 2 Answers

**, see this, this and this. Say, I want to solve \( \Delta \mathbf{u} = \mathbf{f} \), subject to boundary conditions \( \mathbf{u} \cdot \mathbf{t} = g \) and \( \mathbf{u} \cdot \mathbf{n} = 0 \), where \(g\) is some scalar function, e.g., \( g = \nabla\phi \cdot \mathbf{t} \), \( \mathbf{n} \) is the normal vector and \( \mathbf{t} \) is the tangential vector. From the boundary conditions, we have**

*Nitsche's method*\[

\mathbf{u} = (\mathbf{u} \cdot \mathbf{t}) \mathbf{t} = g \mathbf{t} \quad \text{on} \quad \partial \Omega.

\]

Then, we can proceed with the weak formulation as follows:

\[

(\mathbf{f}, \mathbf{v}) = (\Delta \mathbf{u}, \mathbf{v}) = - ( \nabla \mathbf{u}, \nabla \mathbf{v} ) + \int_{\partial \Omega} [\nabla \mathbf{u}] \mathbf{n} \cdot \mathbf{v} \,\mathrm{d} S + \int_{\partial \Omega} [\nabla \mathbf{v}] \mathbf{n} \cdot (\mathbf{u} - g \mathbf{t}) \,\mathrm{d} S,

\]

must hold \( \forall \mathbf{v} \in H^1(\Omega) \), the last term being zero. Following the first link, we add another term which is also zero: \( \beta \int_{\partial \Omega} (\mathbf{u} - g \mathbf{t}) \cdot \mathbf{v} \,\mathrm{d} S \), \( \beta > 0 \) being large enough. The weak formulation then reads:

\[

- ( \nabla \mathbf{u}, \nabla \mathbf{v} ) + \int_{\partial \Omega} [\nabla \mathbf{u}] \mathbf{n} \cdot \mathbf{v} \,\mathrm{d} S + \int_{\partial \Omega} [\nabla \mathbf{v}] \mathbf{n} \cdot \mathbf{u} \,\mathrm{d} S + \beta \int_{\partial \Omega} \mathbf{u} \cdot \mathbf{v} \,\mathrm{d} S =\\ (\mathbf{f}, \mathbf{v}) + \int_{\partial \Omega} [\nabla \mathbf{v}] \mathbf{n} \cdot g \mathbf{t} \,\mathrm{d} S + \beta \int_{\partial \Omega} g \mathbf{t} \cdot \mathbf{v} \,\mathrm{d} S, \quad \forall \mathbf{v} \in H^1(\Omega).

\]

There might be different methods but I tested this approach on a simple domain with

```
from fenics import *
mesh = UnitSquareMesh(24, 24)
V = VectorFunctionSpace(mesh, 'P', 1)
f = Constant((0.0, 0.0))
g = Expression('x[1]', degree=2)
n = FacetNormal(mesh)
t = as_vector([n[1], -n[0]])
u = TrialFunction(V)
v = TestFunction(V)
beta = Constant(10000.0)/mesh.hmax()
a = - inner(grad(u), grad(v))*dx + dot(dot(grad(u), n), v)*ds + dot(u, dot(grad(v), n))*ds + beta*dot(u, v)*ds
L = dot(f, v)*dx + dot(g*t, dot(grad(v), n))*ds + beta*g*dot(t, v)*ds
u = Function(V)
solve(a == L, u)
```

and it seems to be working (even though corners are bit troublesome).

The code works, as you mentioned the corners are troublesome. Any comment for that?

Best, Hamilton

Best, Hamilton

If your mesh is easy, you can write the expression directly. For instance if that condition only occurs at the left side of a unit square, the tangent would be given by (0,-1) and your condition would simply be u[1] = something. In code (only for that side) it would look similar to this:

`bc_left = DirichletBC(V.sub(1), expression_given_by_phi)`

If your mesh is complicated, then a general method for incorporating arbitrary dirichlet boundary conditions exists, called Nitsche's trick/method. You can find it in many places, but I find this reference very good: https://scicomp.stackexchange.com/questions/19910/what-is-the-general-idea-of-nitsches-method-in-numerical-analysis .Best regards!

1, for an easy mesh, for example, a square mesh. expression_given_by_phi can be got easily. What about the complicate geometry, I do not know, how to generate the right term, namely expression_given_by_phi term.Because, expression_given_by_phi = $\left(\nabla\phi\cdot\tau\right)\tau$(∇

`ϕ`·

`τ`)

`τ`. As we know, $\phi$

`ϕ`is defined in domain, and $\tau$

`τ`is defined on boundary. $\phi$

`ϕ`is known by the previous step. How can we calculate $\left(\nabla\phi\cdot\tau\right)\tau$(∇

`ϕ`·

`τ`)

`τ`? I try to export the data of $\nabla\phi$∇

`ϕ`, however, i can not export the normal or tangential directions of boundaries. Could you give me some comment?

2, Your recomment Nitsche's trick/method is valuable for me.

Best, Hamilton

No worries, welcome to the community :).

1 For complicated meshes, the solution should be Nitsche's trick. Another way of doing it would be to create a custom expression for the boundary condition. Something like

```
class CustomExpression(Expression):
def eval(self, x, values):
n = FacetNormal(mesh)(x) # I believe this gives the value at a point
t = as_vector([-n[1],n[0]])
val_phi = nabla_phi(x) # Vector
values = nabla_phi.dot(t)
```

So the idea is that if you are looking to set the values equal to something, you could take the more restricted case of setting \( u \) in the boundary to be exactly \( \nabla \phi \cdot t \), so that they would have the same projection on the tangent. This would work but has some errors. If what you are trying to do is get a field which is the tangent on the boundary, you can solve the Neumann heat equation where the Neumann boundary condition is that the gradient equals the tangent. With that, you get a solution which extends the tangent to the entire domain. With that, you wouldn't need to export any field, because it would be the solution of a preprocessed problem.

2 I love that reference, it is crystal clear

In any case, use the implementation given in the other answer, it looks good and very clean.

Best regs

Thanks for your kind comments. I am fresh one for FEniCS. Some future questions as follows

1, is Custom(Expression) or class(Expression)

2, in my example, $m=\left(\nabla\phi\cdot\tau\right)\cdot\tau$

`m`=(∇

`ϕ`·

`τ`)·

`τ`on the boundaries. where m is a vector and $\phi$

`ϕ`is a scaler. So, the code can be changed into the following one?

#############

Custom forceterm(Expression):

def eval(self, x, values):

n = FacetNormal(mesh)(x) # I believe this gives the value at a point

t = as_vector([-n[1],n[0]])

val_phi = nabla_phi(x) # Vector

values_temp = nabla_phi.dot(t) # scalar

values = values_temp.dot(t)

############

3, if the

*expression_given_by_phi*is defined, how to use this expression, in my sense,

bcs = DirichletBC(Mh, forceterm(Expression),NoSlip)

However, lots of errors accured when the code run. I know my definition is not correct, but i do not know how to make it right .

Best regs

2. If \( \phi \) is a scalar, then \( m \) would be a vector, please recheck your problem

3. I don't understand the question. If you are talking about multiple conditions, refer to this: https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html

Best regards

1, in my problem, $m$

`m`is a vector, the Poisson Equation i expressed above is an example, where $u$

`u`could be a vector. My question is as follows

$m=\left(\nabla\phi\cdot\tau\right)\tau$

`m`=(∇

`ϕ`·

`τ`)

`τ`on the boundaries, where $m$

`m`is a vector and $\phi$

`ϕ`is a scalar. If so, the following code is right?

class CustomExp(Expression):

def eval(self, x, values):

n = FacetNormal(mesh)(x) # I believe this gives the value at a point

t = as_vector([-n[1],n[0]])

val_phi = nabla_phi(x) # Vector

values_temp = nabla_phi.dot(t) # scalar

values = values_temp.dot(t)

############

2, I can understand the tutorial document about the multiple conditions. I do not know how to use CustomExp values.

In my previous code, i just use the following one

bcs = DirichletBC(Mh, CustomExp(phi),NoSlip)

The tutorial tells how to implement the Dirichlet Boundary conditions. For example,

bcs = DirichletBC(Mh, Expression(["0.0","0.0"]),NoSlip)

In my problem, phi is changed every step.

Best, regs.

I try to make it work with the code as follows:

class forceterm(Expression):

def set_phi_value(self, phi):

self.phi = phi

def eval(self, x, values):

n = FacetNormal(mesh)

t = as_vector([-n[1],n[0]])

grad_phi = grad(self.phi) # a vector

values_temp = grad_phi.dot(t) # a scalar

values = values_temp.dot(t) # a vector

and

use the expression I defined with the following code

g = forceterm()

g.set_phi_value = phi_i # phi_i is changed every loop

bcs = DirichletBC(Mh, g,NoSlip)

However, it does not work.

Any comments for that?

Could you help me to give the

*expression_given_by_phi*which represents $\left(\nabla\phi\cdot\tau\right)\tau$(∇

`ϕ`·

`τ`)

`τ`

The code above does not work.

I am eager to know how to implement that.

Best,

Hamilton.

$\Delta u=f$Δ

u=ƒwith the following boundary condition.

$u\cdot\tau=\nabla\phi\cdot\tau$

u·τ=∇ϕ·τ; $u\cdot n=\nabla\phi\cdot n$u·n=∇ϕ·n;where $\tau$

τis the tangential direction of boundaries and $n$nis the normal direction.