How to solve Poisson equation with tangential boundary condition?


279
views
0
4 months ago by
Hi all, 

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
Community: FEniCS Project
Poisson Equation
 $\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$n is the normal direction. 
written 4 months ago by Hamilton  

2 Answers


1
4 months ago by
For simple domains, this shouldn't be a problem. For more complicated ones, the remedy is the Nitsche's method, see thisthis 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
\[
\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).

Thank you so much, you make it clear, and the references let me know more about Nitsche's method.
The code works, as you mentioned the corners are troublesome. Any comment for that?

Best, Hamilton
written 4 months ago by Hamilton  
I think that the main source of the problem are the ill-defined tangential vectors at the corners. It should be ok for domains with a smooth boundary.
written 4 months ago by Adam Janecka  
Appreciate your kind comments. 
Best, Hamilton
written 4 months ago by Hamilton  
One more question, is the weak formulation suitable for Stokes or Navier-Stokes equations?
written 4 months ago by Hamilton  
The Nitsche's method can be used very well for (Navier)-Stokes equations, see the second and third link or https://github.com/MiroK/fenics-nitsche, though the exact weak formulation must be adjusted accordingly.
written 4 months ago by Adam Janecka  
0
4 months ago by
You have two options.

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!
Thanks so much!
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
written 4 months ago by Hamilton  
1

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

written 4 months ago by Nicolás Barnafi  
Hi,
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
written 4 months ago by Hamilton  
1. It is class CustomExp(Expression), I just corrected my answer
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
written 4 months ago by Nicolás Barnafi  
Thanks for your rapid reply.
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.
written 4 months ago by Hamilton  
Hi,
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?
written 4 months ago by Hamilton  
hi, Nicolás Barnafi
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.
written 3 months ago by Hamilton  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »