I’m trying to solve the Poisson equation given as follows (sorry, I don’t how to use Latex here):
delta u = x2 + y2
with the boundary conditions:
u(x,0) = 0
u(0,y) = 0
ux(0,y) = 0
uy(x,0) = 0
Here is my code:
from ngsolve import *
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
source = x**2 + y**2
fes = H1(mesh, order=2, dirichlet='bottom|left')
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes, symmetric=True)
a += -grad(u) * grad(v) * dx
a.Assemble()
f = LinearForm(fes)
f += source * v * dx
f.Assemble()
gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec
Draw(gfu, mesh, "numerical_solution")
From the code above I get an incorrect solution. The exact solution for this equation is (x2y2)/2. What am I doing wrong?
one small Issue: You should consider -Delta u = f. The main issue is, that your solution does not have homogeneous Neumann values (-Grad(u) * n ) on the right and top boundaries, which you imply in your code. You therefore either need to consider inhomogeneous Dirichlet conditions on these boundaries, or implement the inhomogeneous Neumann BC in the forcing term f
Best wishes,
Henry
[code]from ngsolve import *
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
source = x2 + y2
u_ex = - x2 * y2 / 2
fes = H1(mesh, order=2, dirichlet=‘bottom|left’)
u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes, symmetric=True)
a += grad(u) * grad(v) * dx
a.Assemble()
f = LinearForm(fes)
f += source * v * dx
f += - y2 * v * ds(definedon=mesh.Boundaries(‘right’))
f += - x2 * v * ds(definedon=mesh.Boundaries(‘top’))
f.Assemble()