3D Poisson - problem with neumann BC

variational formulation
1

Hi there,

so I want to solve a simple Poisson equation in a unit cube with a neumann BC on one face and a Dirichlet BC on the opposite face.

However the solution looks quite off at the face with the neumann BC. I would expect to have a somewhat uniform distribution.

1
11062×802 57.4 KB

Here is the code:

from fenics import *
import numpy as np

mesh = BoxMesh(Point(-.5,-.5,0), Point(.5,.5,1),20,20,20)

V = FunctionSpace(mesh, 'CG', 1)

def top(x, on_boundary):
    return x[2]>.99  and on_boundary 

def bottom(x, on_boundary):
    return x[2]<.01 and on_boundary

class Inflow(SubDomain):
    def inside(self,x, on_boundary): 
        return bottom(x, on_boundary)

sub_domains = MeshFunction('size_t', mesh, 1)
sub_domains.set_all(0)

inflow_domain = Inflow()
inflow_domain.mark(sub_domains, 1)

ds = Measure('ds', domain=mesh, subdomain_data=sub_domains)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)


a = inner(grad(u),grad(v))*dx
L = Constant(1)*v*ds(1)
 
bc = DirichletBC(V, Constant(1), top)
 
u = Function(V)

solve(a == L, u, bc)

vtkfile = File('solution.pvd')
vtkfile << u
subdomain_file = File("3d/subdomains.pvd")
subdomain_file << sub_domains

Any ideas on why that might happen?

2

On the inflow boundary, you are imposing a Neumann condition, not a Dirichlet condition. If you want to impose a Dirichlet condition weakly, you should consider the method of lagrange multipliers or Nitsches method.

3

What I want is a constant value on one face and a constant inflow on the other. Would that require the method of lagrange multipliers or Nitsches method?

In 2D my implementation works fine.

4

are you trying to impose
u=g or dot(u,n)=g or nu*dot(n, grad(u))=g?

If you are trying to impose something like \int_{\partial\Omega} dot(u,n)ds=g you need to use lagrange multipliers

5

nu*dot(n, grad(u))=g

6

There is a bug in your code:

should be:

sub_domains = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
sub_domains.set_all(0)

as you are marking facets, not edges.

1 Like
7

It works! Thank you very much!