### Poisson Example

329
views
0
10 months ago by
I'm trying solving a simplified version of the poisson example, but it seems I can't get it work.

$u_{xx}+u_{yy}=-xy$uxx+uyy=xy , boundary conditions are  $u\left(x=-1\right)=u\left(x=1\right)=0$u(x=1)=u(x=1)=0.

Exact solution is  $u=\frac{1}{6}xy-\frac{1}{6}x^3y$u=16 xy16 x3y .
I get $u_{\text{max,FEniCS}}\approx0.44$umax,FEniCS0.44 from FEniCS while I expect  $u_{\text{max,exact}}\approx0.67$umax,exact0.67 .

Any idea what could caused this numerical difference? Thanks in advance.

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(2**5, 2**5)

# x and y go from -1 to 1
mesh.coordinates()[:,:] = 2*mesh.coordinates()-1.

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

# Define Dirichlet boundary (x = -1 or x = 1)
def boundary(x):
if near(x[0],1.):
return True
elif near(x[0],-1.):
return True
else:
return False

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('x[0]*x[1]',degree=4)

L = f*v*dx

# Compute solution
u = Function(V,name='FEM')
solve(a == L, u, bc)

# Compute exact solution
u_true = Expression('(1/6.)*x[0]*x[1]-(1/6.)*pow(x[0],3)*x[1]',degree=2)
u_t = interpolate(u_true,V)

# Save solution in VTK format
file = File("poisson.pvd")
file << u

file = File("poisson_t.pvd")
file << u_t

# Plot solution
plot(u, interactive=True)
plot(u_t, interactive=True)​
Community: FEniCS Project

4
10 months ago by
You have not imposed any boundary condition on the top and bottom of the square, y = 1 and y = -1. This means a homogeneous Neumann boundary condition, du/dy = 0 is implied there for your finite element solution. But the function you claim is the exact solution does not satisfy this condition. So your finite element solution is correct, but your "exact solution" is not.
0
10 months ago by
Try to change to V = FunctionSpace(mesh, 'CG', 2) ?
0
10 months ago by
Thank you for your help, Prof. Arnold!