0
$\begingroup$

I'm solving the simple Poisson problem $$-u''(x)=1$$ in the interval $[0,1]$ with $u(0)=u(1)=0$.

I discretised my domain as done here, i.e. with

nodes = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0] and elements = [[0, 1, 2], [2, 3, 4], [4, 5, 6], [6, 7, 8]]

I'm trying to assemble the matrix and the r.h.s by looping over the elements, and I found the same results of here (see pag. 22 for the matrix and load vector)

$$A_K=\frac{1}{3h} \begin{bmatrix} 7 & -8 & 1 \\ -8 & 16 & -8 \\ 1 & -8 & 7 \end{bmatrix} $$

and $$F_K=\frac{h}{6} \begin{bmatrix} 1 \\ 4 \\ 1 \end{bmatrix}$$

This simple codes implements the idea, but my numerical solution is scaled by a factor of $(2h)$: if you set u/(2*h) instead of u, you'll get the right solution exact at grid points).

I think the source of error should be in how I define the step $h$, because those midpoints are confusing me a bit, but I can't figure out how to fix this

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate


a = 0.0
b = 1.0
N = 8
h = (b-a)/(N)

nodes = [0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0]
elements = [[0, 1, 2], [2, 3, 4], [4, 5, 6], [6, 7, 8]]


def assembly_stiff_quadratic():
    A = np.array([[7,-8,1],[-8,16,-8],[1,-8,7]])
    return (A)/(3.0*(h))
    
    
def assembly_load_quadratic():
    return ((h)/6.0)*np.array([1,4,1])



numK = len(elements)
numV = len(nodes)
A = np.zeros([numV,numV])
f = np.zeros(numV)


for K in range (0,numK):
    loc2glb = elements[K]
    AK = assembly_stiff_quadratic()
    print(K,loc2glb)
    bK = assembly_load_quadratic()
    A[np.ix_(loc2glb,loc2glb)] += AK
    f[np.ix_(loc2glb)]+=bK



f = f[1:-1]
A = A[1:-1,1:-1]
u = np.linalg.solve(A,f)


def sol(x):
    return -(0.5)*x*(x-1)

x = np.linspace(0,1,N+1)
plt.plot(x[1:-1],u,'o')
plt.plot(x[1:-1],sol(x[1:-1]),'r')

My question is not about debugging the code, rather about how to do the assembly in this case. I posted the code because I think it is straightforward to read and easy to test.

EDIT:

After the comments, I wrote the following snippet which produces the following order of convergence plot

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate



a = 0.0
b = 1.0
def solvePoisson(N): 
    h = 2*(b-a)/(N)


    elements=[]
    for i in range(0,int(N-1)):

        if i%2==0: elements.append([j for j in range(i,i+3)])
    nodes= np.linspace(a,b,N+1)

    def phi0(x):
        return (1-x)*(1-2*x)

    def phi1(x):
        return 4*x*(1-x)

    def phi2(x):
        return x*(2*x-1)
    
    def rhs(x):
        return np.sin(x)    


    def assembly_stiff_quadratic():
        A = np.array([[7,-8,1],[-8,16,-8],[1,-8,7]])
        return (A)/(3.0*(h))

    numK = len(elements)
    numV = len(nodes)
    A = np.zeros([numV,numV])
    f = np.zeros(numV)



    for K in range (0,numK):
        loc2glb = elements[K]
        AA = nodes[loc2glb[0]]
#        BB = nodes[loc2glb[2]]
        AK = assembly_stiff_quadratic()
        bK = h*np.array([integrate.quadrature(lambda x: rhs(AA+h*x)*phi0(x),0,1)[0],integrate.quadrature(lambda x: rhs(AA+h*x)*phi1(x),0,1)[0],integrate.quadrature(lambda x: rhs(AA+h*x)*phi2(x),0,1)[0]])

        A[np.ix_(loc2glb,loc2glb)] += AK
        f[np.ix_(loc2glb)]+=bK



    f = f[1:-1]
    A = A[1:-1,1:-1]
    u = np.linalg.solve(A,f)
    return u

def sol(x):
    return np.sin(x) - x *np.sin(1)


N = 8
x = np.linspace(0,1,N+1) 
plt.plot(x[1:-1],solvePoisson(N),'o')
plt.plot(x[1:-1],sol(x[1:-1]),'r')

enter image description here

$\endgroup$
12
  • $\begingroup$ Just to make sure: you asked a very similar question and Dr. Bangerth answered it. Why aren't you following his advice? $\endgroup$ – Abdullah Ali Sivas Mar 7 at 14:33
  • $\begingroup$ Actually I think I was following his advice: in each element I have three numbers. Also, my assignment sheet explicitly tells me to use those two arrays nodes and elements @AbdullahAliSivas $\endgroup$ – bobinthebox Mar 7 at 18:13
  • $\begingroup$ What he is suggesting is to have some sort of abstraction layer separating the vertices of the mesh from the nodes of an element. In your case, vertices would be 0, 0.25, 0.5, 0.75 and 1.0 and you would be integrating over [0,0.25], [0.25,05] and so on. The nodes would determine the polynomials you use and they would be sampled more frequently. If you want to keep your code the same, h should be twice the value it is now. $\endgroup$ – Abdullah Ali Sivas Mar 7 at 20:49
  • $\begingroup$ That's the point I'm missing: why should h be doubled? @AbdullahAliSivas Thanks for your time $\endgroup$ – bobinthebox Mar 7 at 21:22
  • $\begingroup$ Think of the mesh as a tessellation of the domain, i.e., $T=\{K\}$. Then on each mesh element $K$, you define a finite element (particularly, Cea's definition, but the modern interpretation is slightly different to allow some exotic elements): the nodes, corresponding polynomials, etc. In your case, 1D, the mesh is [0, 0.25, 0.5, 0.75, 1.0], then on each $K$ you have a finite element with the nodes $[k_0, (k_0+k_1)/2, k_1]$ and the corresponding polynomials $\phi_0,\phi_1,\phi_2$. So the matrix entries you compute are $A^K_{i,j} = \int_K \phi_i\phi_j$ which introduces $h=|K|$. $\endgroup$ – Abdullah Ali Sivas Mar 7 at 21:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.