Poisson Equation FEM

Poisson Equation FEM#

from ngsolve import *
from ngsolve.bem import *
from ngsolve.webgui import Draw
from ngsolve.krylovspace import CG

mesh = Mesh(unit_cube.GenerateMesh(maxh=0.2))
fes = H1(mesh, order=3, dirichlet=".*")

u, v = fes.TnT()
source = 32*(y*(1-y)+x*(1-x))
f = LinearForm(fes)
f += source*v*dx

a = BilinearForm(grad(u)*grad(v)*dx)

a.Assemble()
f.Assemble();

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

Draw (gfu, mesh, clipping={"vec":[0,0,-1], "function": True});

Poisson Equation BEM#

fesL2 = SurfaceL2(mesh, order=2, dual_mapping=False)
u,v = fesL2.TnT()
fesH1 = H1(mesh, order=3)
gfs = GridFunction(fesH1)
gfs.Set(source)

V = LaplaceSL(u*ds) * v*ds
N = LaplaceSL(u*dx)(gfs)

rhs = LinearForm (-N*v.Trace()*ds(bonus_intorder=3)).Assemble()
with TaskManager():
    V = LaplaceSL(u*ds(bonus_intorder=3))*v*ds(bonus_intorder=3)

j = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():  
    CG(mat = V.mat, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, 
       maxsteps=200, initialize=False, printrates=False)
Draw (j);

sol = LaplaceSL(u*ds(bonus_intorder=3))(j) + N
Draw(sol, mesh, clipping={"vec":[0,0,-1], "function": True})

BaseWebGuiScene

Compare solutions#

Draw(sol-gfu, mesh, clipping={"vec":[0,0,-1], "function": True})

BaseWebGuiScene
Contents