Hi,
I’ve been trying to reproduce the 1-D Poisson example from chapter 1 of the Claes Johnson book using Ngsolve with a uniform 1-D mesh of the interval [0,1], a constant forcing function f(x)=1 and an H1 space with order 1.
I was thinking that would be the appropriate space to use, but I’m getting slightly different assembled linear and bilinear forms compared to what I expected. The solution is somewhat different as well.
I’ve been browsing the source code in fem and mylittlengsolve to try to get a better understanding of how the basis functions for a 1-D H1 finite element space are defined, but so far I can’t understand why the solutions don’t match.
I was expecting to get this system and solution (the solution should be unique):
[code]mat
[[ 20. -10. 0. 0. 0. 0. 0. 0. 0. 0.]
[-10. 20. -10. 0. 0. 0. 0. 0. 0. 0.]
[ 0. -10. 20. -10. 0. 0. 0. 0. 0. 0.]
[ 0. 0. -10. 20. -10. 0. 0. 0. 0. 0.]
[ 0. 0. 0. -10. 20. -10. 0. 0. 0. 0.]
[ 0. 0. 0. 0. -10. 20. -10. 0. 0. 0.]
[ 0. 0. 0. 0. 0. -10. 20. -10. 0. 0.]
[ 0. 0. 0. 0. 0. 0. -10. 20. -10. 0.]
[ 0. 0. 0. 0. 0. 0. 0. -10. 20. -10.]
[ 0. 0. 0. 0. 0. 0. 0. 0. -10. 20.]]
rhs
[0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1]
sol
[0.05 0.09 0.12 0.14 0.15 0.15 0.14 0.12 0.09 0.05]
[/code]
but I’m seeing this system. The first and last rows of the matrix and vector are different from what I expected. The solution is a bit flatter:
[code]mat =
Row 0: 0: 10 1: -10
Row 1: 0: -10 1: 20 2: -10
Row 2: 1: -10 2: 20 3: -10
Row 3: 2: -10 3: 20 4: -10
Row 4: 3: -10 4: 20 5: -10
Row 5: 4: -10 5: 20 6: -10
Row 6: 5: -10 6: 20 7: -10
Row 7: 6: -10 7: 20 8: -10
Row 8: 7: -10 8: 20 9: -10
Row 9: 8: -10 9: 20 10: -10
Row 10: 9: -10 10: 10
rhs =
0.05
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.1
0.05
sol =
0
0.045
0.08
0.105
0.12
0.125
0.12
0.105
0.08
0.045
0
[/code]
Here’s the ngsolve code I’m running:
from netgen.meshing import *
# generate a 1D mesh
m = Mesh()
m.dim = 1
nel = 10
pnums = []
for i in range(nel+1):
pnums.append (m.Add (MeshPoint (Pnt(i/nel, 0, 0))))
# add segments
for i in range(nel):
m.Add (Element1D ([pnums[i],pnums[i+1]], index=1))
m.SetMaterial(1,'material')
# add points
m.Add (Element0D (pnums[0], index=1))
m.Add (Element0D (pnums[nel], index=2))
# set boundary condition names
m.SetBCName(0,'left')
m.SetBCName(1,'right')
import ngsolve
from ngsolve import *
ngsmesh = ngsolve.Mesh(m)
print(ngsmesh.GetBoundaries())
# Specify Dirichlet boundary conditions
fes = H1(ngsmesh, order=1, dirichlet='left|right')
#fes = H1(ngsmesh, order=1, dirichlet=[1,2])
print ("freedofs:\n", fes.FreeDofs())
u = fes.TrialFunction() # symbolic object
v = fes.TestFunction() # symbolic object
gfu = GridFunction(fes) # solution
#a = BilinearForm(fes, symmetric=True)
a = BilinearForm(fes)
a += SymbolicBFI(grad(u)*grad(v))
a.Assemble()
print ("mat = \n", a.mat)
# Forcing function
f = CoefficientFunction(1)
lf = LinearForm(fes)
lf += SymbolicLFI(f*v)
lf.Assemble()
print ("rhs = \n", lf.vec)
u = GridFunction(fes)
u.vec.data = a.mat.Inverse(fes.FreeDofs()) * lf.vec
print ("sol =\n", u.vec)
pnts = [i/100 for i in range(101)]
pnts_vals = [ (x,u(x)) for x in pnts if ngsmesh.Contains(x)]
%matplotlib inline
import matplotlib.pyplot as plt
pnts,vals = zip(*pnts_vals)
plt.plot(pnts,vals, "-*")
plt.show()
This code generates what I expected.
from numpy import eye, diag, ones
from scipy.linalg import solve
h = 0.1
A = ( 2*eye(10)-diag(ones(9),1)-diag(ones(9),-1) )/h
b = h*ones(10)
u = solve(A,b)
Thanks for any insight on this!
Best,
Dow
