Poisson Equation

\[\begin{split}\begin{cases} -\nabla^2 u = f(x) &\text{ in } \Omega \\ u = 0 & \text{ on } \partial\Omega \end{cases}\end{split}\]

1 Dimensional Case

Example
>>> import numpy as np
>>> from mozart.mesh.rectangle import interval
>>> from mozart.poisson.fem.interval import solve
>>> f = lambda x: np.pi ** 2 * np.sin(np.pi * x)
>>> u_D = lambda x: np.zeros_like(x)
>>> nrElems, degree = (7, 1)
>>> c4n, n4e, n4db, ind4e = interval(0, 1, nrElems, degree)
>>> u = solve(c4n, n4e, n4db, ind4e, f, u_D, degree)
>>> u
array([ 0.        ,  0.42667492,  0.76884164,  0.95872984,  0.95872984,
0.76884164,  0.42667492,  0.        ])
Poisson 1D Figure
mozart.poisson.fem.interval.DVandermondeM1D(degree, r)[source]

Initialize the derivative of modal basis (i) at (r) at order degree

Paramters
  • degree (int32) : Polynomial degree
  • r (float64 array) : points
Returns
  • DVr (float64 array) : Differentiate Vandermonde matrix
Example
>>> N = 3
>>> from mozart.poisson.fem.interval import VandermondeM1D
>>> r = np.linspace(-1,1,N+1)
>>> DVr = DVandermondeM1D(N,r)
>>> print(DVr)
array([[  0.        ,   1.22474487,  -4.74341649,  11.22497216],
[  0.        ,   1.22474487,  -1.58113883,  -1.24721913],
[  0.        ,   1.22474487,   1.58113883,  -1.24721913],
[  0.        ,   1.22474487,   4.74341649,  11.22497216]])
mozart.poisson.fem.interval.Dmatrix1D(degree, r, V)[source]

Initialize the derivative of modal basis (i) at (r) at order degree

Paramters
  • degree (int32) : Polynomial degree
  • r (float64 array) : points
Returns
  • DVr (float64 array) : Differentiate Vandermonde matrix
Example
>>> N = 3
>>> from mozart.poisson.fem.interval import Dmatrix1D
>>> r = np.linspace(-1,1,N+1)
>>> Dr = Dmatirx1D(N,r)
>>> print(Dr)
array([[-2.75,  4.5 , -2.25,  0.5 ],
[-0.5 , -0.75,  1.5 , -0.25],
[ 0.25, -1.5 ,  0.75,  0.5 ],
[-0.5 ,  2.25, -4.5 ,  2.75]])
mozart.poisson.fem.interval.VandermondeM1D(degree, r)[source]

Initialize the 1D Vandermonde matrix, \(V_{i,j} = \phi_j(r_i)\)

Paramters
  • degree (int32) : Polynomial degree
  • r (float64 array) : points
Returns
  • V1D (float64 array) : 1D Vandermonde matrix
Example
>>> N = 3
>>> from mozart.poisson.fem.interval import VandermondeM1D
>>> r = np.linspace(-1,1,N+1)
>>> V1D = VandermondeM1D(N,r)
>>> print(V1D)
array([[ 0.70710678, -1.22474487,  1.58113883, -1.87082869],
[ 0.70710678, -0.40824829, -0.52704628,  0.76218947],
[ 0.70710678,  0.40824829, -0.52704628, -0.76218947],
[ 0.70710678,  1.22474487,  1.58113883,  1.87082869]])
mozart.poisson.fem.interval.computeError(c4n, n4e, ind4e, exact_u, exact_ux, approx_u, degree, degree_i)[source]

Computes L^2-error and semi H^1-error between exact solution and approximate solution.

Parameters
  • c4n (float64 array) : coordinates for nodes
  • n4e (int32 array) : nodes for elements
  • ind4e (int32 array) : indices for elements
  • exact_u (lambda) : exact solution
  • exact_ux (lambda) : derivative of exact solution
  • approx_u (float64 array) : approximate solution
  • degree (int32) : Polynomial degree
  • degree_i (int32) : Polynomial degree for interpolation
Returns
  • L2error (float64) : L^2 error between exact solution and approximate solution.
  • sH1error (float64) : semi H^1 error between exact solution and approximate solution.
Example
>>> N = 2
>>> from mozart.mesh.rectangle import interval 
>>> c4n, n4e, n4db, ind4e = interval(0, 1, 4, 2)
>>> f = lambda x: np.pi ** 2 * np.sin(np.pi * x)
>>> u_D = lambda x: np.zeros_like(x)
>>> from mozart.poisson.fem.interval import solve_p
>>> x = solve_p(c4n, n4e, n4db, ind4e, f, u_D, N)
>>> from mozart.poisson.fem.interval import computeError
>>> exact_u = lambda x: np.sin(np.pi * x)
>>> exact_ux = lambda x: np.pi * np.cos(np.pi * x)
>>> L2error, sH1error = computeError(c4n, n4e, ind4e, exact_u, exact_ux, x, N, N+3)
>>> L2error
0.0020225729623142077
>>> sH1error
0.05062779815975444
mozart.poisson.fem.interval.getMatrix(degree)[source]

Get FEM matrices on the reference domain I = [-1, 1]

Paramters
  • degree (int32) : degree of polynomial
Returns
  • M_R (float64 array) : Mass matrix on the reference domain
  • S_R (float64 array) : Stiffness matrix on the reference domain
  • D_R (float64 array) : Differentiation matrix on the reference domain
mozart.poisson.fem.interval.solve(c4n, n4e, n4db, ind4e, f, u_D, degree)[source]

Computes the coordinates of nodes and elements.

Parameters
  • c4n (float64 array) : coordinates for nodes
  • n4e (int32 array) : nodes for elements
  • n4db (int32 array) : nodes for Dirichlet boundary
  • ind4e (int32 array) : indices for elements
  • f (lambda) : source term
  • u_D (lambda) : Dirichlet boundary condition
  • degree (int32) : Polynomial degree
Returns
  • x (float64 array) : solution
Example
>>> N = 2
>>> from mozart.mesh.rectangle import interval 
>>> c4n, n4e, n4db, ind4e = interval(0, 1, 4, 2)
>>> f = lambda x: np.ones_like(x)
>>> u_D = lambda x: np.zeros_like(x)
>>> from mozart.poisson.fem.interval import solve
>>> x = solve(c4n, n4e, n4db, ind4e, f, u_D, N)
>>> x
array([ 0.       ,  0.0546875,  0.09375  ,  0.1171875,  0.125    ,
   0.1171875,  0.09375  ,  0.0546875,  0.       ])