Forum Message

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

1-D Poisson with constant forcing function

ddrake
Topic Author
Offline
Premium Member

5 years 10 months ago - 5 years 10 months ago #554 by ddrake

1-D Poisson with constant forcing function was created by ddrake

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]

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

Here's the ngsolve code I'm running:

Code:

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.

Code:

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

Last edit: 5 years 10 months ago by ddrake. Reason: correct error in terminology

The following user(s) said Thank You: liubenyuan

joachim
Offline
Administrator

5 years 10 months ago #556 by joachim

Replied by joachim on topic 1-D Poisson with constant forcing function

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

The following user(s) said Thank You: ddrake

ddrake
Topic Author
Offline
Premium Member

5 years 10 months ago #557 by ddrake

Replied by ddrake on topic 1-D Poisson with constant forcing function

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 :-$

Time to create page: 0.113 seconds

Login

NGSolve on Twitter