1-D Poisson with constant forcing function

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))))

for i in range(nel):

m.SetMaterial(1,'material')

# 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.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

Hi Dow,

I assume you want to have 10 elements, i.e. 9 interior points and 11 points in total.
NGSolve computs the matrix and right hand side for the Neumann problem (11 points). The boundary constraints are taken into account by the Inverse (FreeDofs()).

Your matrix a la Claes Johnson is 10 x 10, but should be 9 x 9 (number of interior points).

btw, the exact solution in x=0.5 is 1/8 = 0.125

Joachim

Thank you very much, Dr. Schoeberl!

I need to teach a 2 hour mini course on FEM and Ngsolve/Netgen next week. It will work better if I know what I’m talking about