Note: This tutorial was generated from an IPython notebook that can be downloaded here.

2D Poisson Equation with Triangular Elements

This tutorial was made with the following version of FEMpy:

import FEMpy
FEMpy.__version__

'1.0'

FEMpy can solve 2D domains using similar inputs as in the 1D case. Here we will solve the following problem

\[-\nabla \cdot \left(c(x,y) \nabla u(x,y)\right) = f(x,y); (x,y) \in \Omega = [-1, 1] \times [-1, 1]\]

where \(c(x,y) = 1\) and \(f(x,y) = -2 e^{x + y}\) and with boundary conditions

\[u(x,y) = e^{-1 + y} \text{ for } x = -1, y \in \delta\Omega \setminus \Gamma_1\]

\[u(x,y) = e^{1 + y} \text{ for } x = 1, y \in \delta\Omega \setminus \Gamma_1\]

\[u(x,y) = e^{x + 1} \text{ for } y = 1, x \in \delta\Omega \setminus \Gamma_1\]

\[\nabla u(x,y) \cdot \hat{\mathbf{n}} = -e^{x - 1} \text{ for } y = -1, x \in \Gamma_1 \subseteq \delta\Omega\]

Let us define our necessary functions

import numpy as np


def coefficient_function(coords):
    return 1


def source_function(coord):
    x, y = coord
    return -2 * np.exp(x + y)


def dirichlet_function(coord):
    x, y = coord
    if x == -1:
        return np.exp(-1 + y)
    elif x == 1:
        return np.exp(1 + y)
    elif y == 1:
        return np.exp(x + 1)
    elif y == -1:
        return np.exp(x - 1)


def neumann_function(coord):
    x, y = coord
    if y == -1:
        return -np.exp(x - 1)

Let’s choose a linear basis for our problem and set up our grid to have step sizes of \(h_1 = h_2 = 1/4.\)

basis = FEMpy.TriangularBasis2D('linear')
mesh = FEMpy.TriangularMesh2D(-1, 1, -1, 1, h1=1/4, h2=1/4, basis_type='linear')

Before we can set up our boundary conditions, we will need to define our boundary node and edge types. These can most easily be done by a list of strings indicating the boundary condition type for each boundary node and edge, respectively.

boundary_node_types = ['dirichlet', *['neumann']*7, *['dirichlet']*9, *['dirichlet']*8, *['dirichlet']*7]
boundary_edge_types = [*['neumann']*8, *['dirichlet']*8, *['dirichlet']*8, *['dirichlet']*8]

Now, we can define our boundary conditions, remembering to include the test basis function and the coefficient function since we have Neumann boundary conditions.

bcs = FEMpy.BoundaryConditions2D(mesh, boundary_node_types, boundary_edge_types,
                                 test_basis=basis,
                                 dirichlet_fun=dirichlet_function,
                                 neumann_fun=neumann_function,
                                 coeff_fun=coefficient_function)

We can input our mesh, basis functions, and boundary conditions into our Poisson equation and call our solver:

poisson_eq = FEMpy.Poisson2D(mesh, fe_test_basis=basis, fe_trial_basis=basis, boundary_conditions=bcs)
poisson_eq.solve(coefficient_function, source_function)

array([0.13533528, 0.17377394, 0.22313016, 0.2865048 , 0.36787944,
       0.47236655, 0.60653066, 0.77880078, 1.        , 0.17325304,
       0.22277292, 0.28625562, 0.36770659, 0.47224895, 0.60645367,
       0.77875452, 0.99997833, 1.28402542, 0.22221544, 0.28584593,
       0.36741281, 0.47204193, 0.60631008, 0.77865671, 0.99991359,
       1.28398499, 1.64872127, 0.28526803, 0.36698257, 0.4717327 ,
       0.60609224, 0.77850475, 0.99980769, 1.28391054, 1.64866766,
       2.11700002, 0.36639009, 0.47130853, 0.60579719, 0.7783016 ,
       0.99966716, 1.28381118, 1.64859408, 2.11694086, 2.71828183,
       0.47072395, 0.60541827, 0.77805796, 0.99950801, 1.28370351,
       1.64851667, 2.11687951, 2.71822603, 3.49034296, 0.60490964,
       0.77779468, 0.99936654, 1.28362138, 1.64846308, 2.11683828,
       2.71818738, 3.49029942, 4.48168907, 0.77755761, 0.99934246,
       1.28364454, 1.64848877, 2.11685491, 2.71819212, 3.49029095,
       4.48166518, 5.75460268, 1.        , 1.28402542, 1.64872127,
       2.11700002, 2.71828183, 3.49034296, 4.48168907, 5.75460268,
       7.3890561 ])

We can also look at the \(L^\infty\) and \(L^2\) norm errors as well as the \(H^1\) semi-norm error associated with our solution as compared against the analytical solution of

\[u(x,y) = e^{x + y}\]

def analytic_solution(coord):
    x, y = coord
    return np.exp(x + y)


def dx_analytic_solution(coord):
    x, y = coord
    return np.exp(x + y)


def dy_analytic_solution(coord):
    x, y = coord
    return np.exp(x + y)


L_inf_err = poisson_eq.l_inf_error(analytic_solution)
L_2_err = poisson_eq.l2_error(analytic_solution)
H_1_err = poisson_eq.h1_seminorm_error((dx_analytic_solution, dy_analytic_solution))

from IPython.display import display, Math
display(Math('\|L^\infty\| = {l_inf:.3e}'.format(l_inf=L_inf_err)))
display(Math('\|L^2\| = {l_2:.3e}'.format(l_2=L_2_err)))
display(Math('|H^1| = {h_1:.3e}'.format(h_1=H_1_err)))

$$\|L^\infty\| = 3.176e+00$$

$$\|L^2\| = 7.357e-04$$

$$|H^1| = 1.417e-02$$