normalization (average zero) condition for Neumann problem did not work

Asked by Daniel Bare on 2013-01-30

Hello,

 I have a simple (solvable) Poisson-Neumann problem. I solved it and made an a posteriori 'normalization', by subtracting the average of the solution. But the function still has no zero average.

Below the code, thanks!!!

from dolfin import *

# Create mesh and define function space
mesh = UnitSquare(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

n=FacetNormal(mesh)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("-2.0*x[0]")
g = Expression(("-0.33*0.5*(pow(x[0],2)+pow(x[1],2))","-0.33*x[0]*x[1]"))
h = Expression("x[0]")
a = inner(grad(u), grad(v))*dx
L = f*v*dx+dot(g,n)*v*ds

# Compute solution
u = Function(V)
solve(a == L, u)

print assemble(u * dx, mesh=mesh)

u_array = u.vector().array() - assemble(u * dx, mesh=mesh)
u.vector()[:] = u_array

average_u = assemble( u * dx, mesh=mesh)

print average_u

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-01-30
Last reply:
2013-01-30
Anders Logg (logg) said : #1

Try

u_array = u.vector().array() - assemble(u * dx, mesh=mesh) / assemble(Constant(1.0)*dx, mesh=mesh)

You might also be interested in looking at two other options:

1. Using Lagrange multipliers. That way you will also get a
nonsingular system to solve. Look at the demo

demo/undocumented/neumann-poisson/python/demo_neumann-poisson.py

2. The function normalize().

--
Anders

On Wed, Jan 30, 2013 at 12:41:08AM -0000, Daniel Bare wrote:
> New question #220505 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/220505
>
>
> Hello,
>
>
>
> I have a simple (solvable) Poisson-Neumann problem. I solved it and made an a posteriori 'normalization', by subtracting the average of the solution. But the function still has no zero average.
>
> Below the code, thanks!!!
>
> from dolfin import *
>
> # Create mesh and define function space
> mesh = UnitSquare(32, 32)
> V = FunctionSpace(mesh, "Lagrange", 1)
>
> n=FacetNormal(mesh)
>
> # Define variational problem
> u = TrialFunction(V)
> v = TestFunction(V)
> f = Expression("-2.0*x[0]")
> g = Expression(("-0.33*0.5*(pow(x[0],2)+pow(x[1],2))","-0.33*x[0]*x[1]"))
> h = Expression("x[0]")
> a = inner(grad(u), grad(v))*dx
> L = f*v*dx+dot(g,n)*v*ds
>
> # Compute solution
> u = Function(V)
> solve(a == L, u)
>
> print assemble(u * dx, mesh=mesh)
>
> u_array = u.vector().array() - assemble(u * dx, mesh=mesh)
> u.vector()[:] = u_array
>
>
> average_u = assemble( u * dx, mesh=mesh)
>
> print average_u
>

Can you help with this problem?

Provide an answer of your own, or ask Daniel Bare for more information if necessary.

To post a message you must log in.