## setup and notebook configuration


```python
# scientific python stack
import numpy as np
import scipy as sp
import sympy as sym
import orthopy, quadpy
```


```python
# matplotlib, plotting setup
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.tri as mtri               # delaunay triangulation
from mpl_toolkits.mplot3d import Axes3D     # surface plotting
import seaborn as sns                       # nice plotting defaults
import cmocean as cmo                       # ocean colormaps
sym.init_printing(use_latex='mathjax')
sns.set()
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
```

    /Users/corbin/virtual_envs/FEMexperiment_env/lib/python3.6/site-packages/cmocean/tools.py:76: MatplotlibDeprecationWarning: The is_string_like function was deprecated in version 2.1.
      if not mpl.cbook.is_string_like(rgbin[0]):



```python
%load_ext autoreload
%autoreload 2
```


```python
# local imports
import src.fem_base.master.mk_basis_nodal as mbn
import src.fem_base.master.mk_master as mkm
import src.fem_maps.fem_map as fem_map
```

# creation of 1D nodal bases

We define the master 1D element as $I\in[-1, 1]$


```python
b = mbn.Basis_nodal(order=1, dim=1, element=0)
b.plot_elm()
```


    

    



```python
fig, ax = plt.subplots(1, 5, figsize = (10, 1))
xx = np.linspace(-1, 1, 100)
pts = [[x, 0, 0] for x in xx]

for order in range(5):
    b = mbn.Basis_nodal(order=order, dim=1, element=0)
    yy = b.eval_at_pts(pts=pts)
    for basis_fn in range(b.nb):
        ax[order].plot(xx, yy[:, basis_fn])
    ax[order].set_title(r'$p = {}$'.format(order))
```


    

    


## construction of vandermonde matrices


```python
x = np.linspace(-1, 1, 100)
n_polys = 4
vals = orthopy.line_segment.tree_jacobi(x, n=n_polys, alpha=0, beta=0, standardization='normal')
for i in range(n_polys):
    plt.plot(x, vals[i], label='P_{}'.format(i))
plt.legend()

plt.title('normalized Legendre polynomials')
plt.show()
```


    

    


These polynomials agree with the explicitly listed polynomials in Hesthaven, so we know that they are orthonormalized correctly.


```python
def Jacobi_Poly(r, alpha, beta, N):
    """ wraps orthopy to return Jacobi polynomial """
    return orthopy.line_segment.tree_jacobi(r, n=N-1, alpha=alpha, beta=beta, standardization='normal')

def P_tilde(r, N):
    P = np.zeros((len(r), N))
    polyvals = Jacobi_Poly(r, alpha=0, beta=0, N=N)
    for j in range(N):
        P[:, j] = polyvals[j]
    return P.T
    
def Vandermonde1D(N, x):
    """ initialize 1D vandermonde Matrix Vij = phi_j(x_i)"""
    V1D = np.zeros((len(x), N))
    JacobiP = Jacobi_Poly(x, alpha=0, beta=0, N=N)
    for j, polyvals in enumerate(JacobiP):
        V1D[:, j] = polyvals
    return V1D
```


```python
def LegendreGaussLobatto(N):
    GL = quadpy.line_segment.GaussLobatto(N, a=0., b=0.)
    return GL.points, GL.weights

def GaussLegendre(N):
    GL = quadpy.line_segment.GaussLegendre(N)
    return GL.points, GL.weights
```

An important conceptual point is that the Vandermonde matrix here is NOT the shape function matrix, it's the Vandermonde matrix of the Orthonormal polynomial basis. We will see this later as we have to create the shape function matrices.

## properties / conditioning of vandermonde matrices


```python
equi_det, LGL_det = [], []
for N in range(2, 35):
    nb = N + 1
    equi_pts = np.linspace(-1, 1, nb)
    V = Vandermonde1D(nb, equi_pts)
    equi_det.append(np.linalg.det(V))
    
    LGL_pts, _ = LegendreGaussLobatto(nb)
    V = Vandermonde1D(nb, LGL_pts)
    LGL_det.append(np.linalg.det(V))
```


```python
plt.semilogy(list(range(2, 35)), equi_det, label='equidistant')
plt.semilogy(list(range(2, 35)), LGL_det, label='LGL nodes')
plt.legend()
plt.show()
```


    

    


This result agrees with Hesthaven.


```python
# construct generic lagrange interpolant
from scipy.interpolate import lagrange

def lagrange_polys(pts):
    lagrange_polys = []
    for i, pt in enumerate(pts):
        data = np.zeros_like(pts)
        data[i] = 1
        lagrange_polys.append(lagrange(pts, data))
    return lagrange_polys

def lagrange_basis_at_pts(lagrange_polys, eval_pts):
    """ evaluates lagrange polynomials at eval_pts"""
    result = np.zeros((len(lagrange_polys) ,len(eval_pts)))
    for i, poly in enumerate(lagrange_polys):
        result[i, :] = lagrange_polys[i](eval_pts)
    return result
```

plot lagrange polys over equally spaced vs LGL points


```python
N = 5
lp = np.linspace(-1, 1, N)
lpolys = lagrange_polys(lp)

vN = 100
view_pts = np.linspace(-1, 1, vN)
li = lagrange_basis_at_pts(lpolys, view_pts)
plt.plot(view_pts, li.T)
plt.title('lagrange polynomials over equally spaced points')
plt.show()
```


    

    



```python
N = 5
lp, _ = LegendreGaussLobatto(N)
lpolys = lagrange_polys(lp)

vN = 100
view_pts = np.linspace(-1, 1, vN)
li = lagrange_basis_at_pts(lpolys, view_pts)
plt.plot(view_pts, li.T)
plt.title('lagrange polynomials over LGL points')
plt.show()
```


    

    


So beautiful! By moving the lagrange data points to the nodal points, our basis functions don't exceed 1, unlike in the above plot, where we are already seeing a slight Runge phenomenon.

## vandermonde relations

### relationship between vandermonde $V$, basis polynomials $\tilde{\mathbf{P}}$, and lagrange basis functions (shape functions) $\ell$

Hesthaven makes the claim that $V^T \mathbf{\ell}(r) = \tilde{\mathbf{P}}(r)$ in (3.3).

In Hesthaven's notation, $N$ denotes the polynomial order, $N_p$ denotes the number of nodal points (we would call $nb$), and let's call the number of "view points" `xx`, which are arbitrary.

Then the shapes of the Hesthaven structures are:
- $\mathbf{\ell}$, $\tilde{\mathbf{P}}$, $V$ are all (`nb`, `xx`)
- $V^T \ell$ is (`xx`, `nb`) x (`nb`, `xx`) $\rightarrow$ (`xx`, `xx`) where rows contain the values of polynomials $\tilde{\mathbf{P}}$

This works for either equidistant points or the LGL points.


```python
N = 5
lp, _ = LegendreGaussLobatto(N) 
#lp = np.linspace(-1, 1, N)
view_pts = np.linspace(-1, 1, 50)

l_polys = lagrange_polys(pts=lp)
ℓ = lagrange_basis_at_pts(l_polys, eval_pts=view_pts)
V = Vandermonde1D(N=len(view_pts), x=lp)
P = np.dot(V.T, ℓ)

# plot the result
plt.plot(view_pts, ℓ.T, '--')
plt.plot(view_pts, P[0:3,:].T)
plt.show()
```


    

    


We see that indeed we recover the Legendre polynomials.

More directly, we can invert the relation to find that

$$\ell = (V^T)^{-1} \tilde{\mathbf{P}}$$

which allows us to create our nodal shape functions.


```python
nb = 4
nodal_pts, _ = LegendreGaussLobatto(nb)
view_pts = np.linspace(-1, 1, 50)

# create the Vandermonde, P matrices
V = Vandermonde1D(N=nb, x=nodal_pts)
Vti = np.linalg.inv(V.T)
P = P_tilde(r=view_pts, N=nb)
print('shape of Vandermonde: {}'.format(V.shape))
print('shape of P: {}'.format(P.shape))

yy = np.dot(Vti, P)
plt.plot(view_pts, yy.T)
plt.title('nodal shape functions generated from orthogonal basis polynomials')
plt.show()
```

    shape of Vandermonde: (4, 4)
    shape of P: (4, 50)



    

    


### relationship between vandermonde $V$ and mass matrix

We can build on the relationship developed in the section above to form the mass matrix for a nodal basis. We note that 

$M_{ij} = \int_{-1}^{1}\ell_i(r)\, \ell_j(r) \,dr = (\ell_i, \ell_j)_I$, and if we expand out $\ell = (V^T)^{-1}\tilde{\mathbf{P}}$, it turns out (page 51)

$$M = (V V^T)^{-1}$$

because of the orthogonal nature of our choice of basis function; the implication is that we can compute the integrals over the master element without the explicit need for quadrature points or weights. Note first that $\phi_i(\xi) = \sum_{n=1}^{nb} (V^T)_{in}^{-1} \tilde{P}_{n-1}(\xi)$. Then

\begin{align}
M_{ij} &= \int^{1}_{-1} \phi_i(\xi)\,\phi_j(\xi)\,d\xi 
      = \int^{1}_{-1}\left[\sum_{k=1}^{nb} (V^T)_{ik}^{-1} \tilde{P}_{k-1}(\xi)
      \sum_{m=1}^{nb} (V^T)_{jm}^{-1} \tilde{P}_{m-1}(\xi) \right]\, d\xi \\
      &= \sum_{k=1}^{nb} \sum_{m=1}^{nb} (V^T)_{ik}^{-1} 
       (V^T)_{jm}^{-1} \int^{1}_{-1}\tilde{P}_{k-1}(\xi) \tilde{P}_{m-1}(\xi)
       =\sum_{k=1}^{nb} \sum_{m=1}^{nb} (V^T)_{ik}^{-1}
       (V^T)_{jm}^{-1} \delta_{km} \\
       &=\sum_{k=1}^{nb} (V^T)_{im}^{-1} 
       (V^T)_{jm}^{-1} = \sum_{k=1}^{nb} (V^T)_{mi}^{-1} (V)_{mj}^{-1} \\
       &= (V^{T})^{-1} V^{-1} = (VV^T)^{-1}
\end{align}

Where note we've used the cute trick that $\int_{-1}^1 \tilde{P}_m \tilde{P}_n = \delta_{mn}$, since we chose an __orthonormal__ modal basis. Orthogonal wouldn't have done it, but an orthonormal modal basis has this property.

We can check this relation against the more traditional way of constructing the mass matrix with quadrature. `master.shap` has dimensions of (`n_quad`, `nb`)


```python
order = 3
m1d = mkm.Master_nodal(order=order, dim=1, element=0)
xq, wq = m1d.cube_pts, m1d.cube_wghts
shap = m1d.shap
shapw = np.dot(np.diag(wq), m1d.shap)
M_quadrature = np.dot(shap.T, shapw)


Np = order + 1
nodal_points, _ = LegendreGaussLobatto(Np)
V = Vandermonde1D(N=Np, x=nodal_points)
M_vand = np.linalg.inv(np.dot(V, V.T))

# this will throw an error if not correct
assert(np.allclose(M_quadrature, M_vand))
```

## efficient computation of derivatives of the basis functions

### derivatives of Legendre polynomials

In order to compute the derivatives of the shape functions (which are expressed via the vandermonde matrix $V$), we must take the derivatives with respect to the orthogonal basis polynomials. There is an identity (Hesthaven, p. 52)

$$ \frac{d \tilde{P}_n}{d r} = \sqrt{n(n+1)}\,\tilde{P}^{(1,1)}_{n-1}$$

This is in contrast to directly differentiating either the coefficients of $\tilde{P}$ or more directly the nodal shape functions $\ell$ if the explicit polynomial form is known (like in `scipy`, but this becomes trickier in multiple dimensions). As it turns out, the first approach is a very efficient way to compute these operators.


```python
def Jacobi_Poly_Derivative(r, alpha, beta, N):
    """ take a derivative of Jacobi Poly, more general than above
    copy the format of orthopy (list of arrays)
    """
    dp = [np.zeros_like(r)]
    Jacobi_P = Jacobi_Poly(r, alpha + 1, beta + 1, N)
    for n in range(1, N+1):
        gamma = np.sqrt(n * (n + alpha + beta + 1))
        dp.append(gamma * Jacobi_P[n-1])
    return dp

#def dP_tilde(r, N):
#    P = np.zeros((len(r), N))
#    polyvals = Jacobi_Poly_Derivative(r, alpha=0, beta=0, N=N)
#    for j in range(N):
#        P[:, j] = polyvals[j]
#    return P
```

We can examine some of the derivatives of the Legendre polynomials.


```python
# some unit testing 
# first jacobi poly is const, so der should be 0
xx = np.linspace(-1, 1, 50)
jpd = Jacobi_Poly_Derivative(xx, alpha=0, beta=0, N=3)
for i, polyder in enumerate(jpd):
    plt.plot(xx, polyder, label=r'$P_{}^\prime(x)$'.format(i))
plt.legend()
plt.show()
```


    

    


These look good. The derivative of the first Legendre polynomial is analytically 0, $P_1^\prime = \sqrt{3/2}$, $P_2^\prime$ should be linear, $P_3^\prime$ should be quadratic.

### discrete derivative operators

We can declare the derivative Vandermonde matrix, and invert it in the same manner to obtain the derivatives of the nodal shape functions.

This works because 
$$V^T \ell = P \Rightarrow V^T \frac{d}{dx} \ell = \frac{d}{dx}P$$

Hence $$V_r\equiv V^T D_r^T, \qquad {V_r}_{(ij)} = \frac{d \tilde{P}_j(r_i)}{d x} $$

and finally $D_r = V_r V^{-1}$ (see Hesthaven, p. 53), as well as $S = M D_r $, where $S_{ij} = \left(\phi_i, \frac{d\phi_j}{dx}\right)_I$, and where $M$ is the mass matrix.


```python
def GradVandermonde1D(N, x):
    Vr = np.zeros((len(x), N))
    dJacobi_P = Jacobi_Poly_Derivative(x, alpha=0, beta=0, N=N-1)
    for j, polyder in enumerate(dJacobi_P):
        Vr[:,j] = polyder
    return Vr
```


```python
p = 3
nb = p+1
nodal_pts, _ = LegendreGaussLobatto(nb)
#nodal_pts = np.linspace(-1, 1, nb)
view_pts = np.linspace(-1, 1, 50)

# grad vandermonde
V  = Vandermonde1D(N=nb, x=nodal_pts)
Vr = GradVandermonde1D(N=nb, x=view_pts)
Vi = np.linalg.inv(V)
Dr = np.dot(Vr, Vi)

print('shape Vr:    {}'.format(Vr.shape))
print('shape V inv: {}'.format(Vi.shape))
print('shape Dr:    {}'.format(Dr.shape))

# shape functions
V = Vandermonde1D(N=nb, x=nodal_pts)
Vti = np.linalg.inv(V.T)
P = P_tilde(r=view_pts, N=nb)
shap = np.dot(Vti, P)

# shape functions at view points
fig, ax = plt.subplots(1, 2, figsize=(8, 4))
ax[0].plot(view_pts, shap.T, '--')
ax[0].set_title(r'nodal shape functions $\phi$')
ax[1].plot(view_pts, Dr)
ax[1].set_title(r'derivatives of nodal shape functions $\frac{d \phi}{dx}$')
plt.show()
```

    shape Vr:    (50, 4)
    shape V inv: (4, 4)
    shape Dr:    (50, 4)



    

    


As a remark, we can once again show the effect of using Legendre Gauss Lobatto points vs equally spaced nodal points.


```python
N = 8
nb = N+1
nodal_pts_LGL, _ = LegendreGaussLobatto(nb)
nodal_pts_equi = np.linspace(-1, 1, nb)
view_pts = np.linspace(-1, 1, 100)

# shape functions at view points
fig, ax = plt.subplots(figsize=(8, 8))
ax.set_yticks([])
ax.set_xticks([])
labels = ['LGL nodes', 'uniform nodes']
for idx, nodal_pts in enumerate([nodal_pts_LGL, nodal_pts_equi]):
    # grad vandermonde
    V  = Vandermonde1D(N=nb, x=nodal_pts)
    Vr = GradVandermonde1D(N=nb, x=view_pts)
    Vi = np.linalg.inv(V)
    Dr = np.dot(Vr, Vi)

    # shape functions
    V = Vandermonde1D(N=nb, x=nodal_pts)
    Vti = np.linalg.inv(V.T)
    P = P_tilde(r=view_pts, N=nb)
    shap = np.dot(Vti, P)
    
    # plot
    ax = fig.add_subplot(2, 2, idx*2+1)
    ax.plot(view_pts, shap.T)
    ax.set_yticks([0, 1])
    ax.set_title(r' $\phi$, {}'.format(labels[idx]))
    ax = fig.add_subplot(2, 2, idx*2+2)
    ax.plot(view_pts, Dr)
    ax.set_title(r'$\f{}$, {}'.format('rac{d \phi}{dx}',labels[idx]))
plt.subplots_adjust(wspace=0.2, hspace=0.2)
fig.suptitle(r'$\phi$ and $\frac{d\phi}{d x}$, LGL vs uniformly-spaced nodes')
plt.show()
```


    

    


### remarks on discrete derivative operators

Suppose we compute the derivative matrix $D_r$ at the nodal points for some order $p$. There are some interesting properties to understand about these derivative matrices.

#### annihilation of constant vectors

Note that if we represent a function nodally, i.e., $u = c_i \phi_i(x)$, then

$$ \frac{du}{dx} = \frac{d}{dx}(c_i \phi_i(x)) = c_i \frac{d\phi_i}{dx} $$

Therefore, if we want to discretely take a derivative of a function (we'll represent the function on the master element for now, but trivially, we could map it to some other region), it suffices to multiply the derivative operator with the nodal vector $D_r u$. It should be clear, then, that $D_r$ will annihilate any constant vector.


```python
p = 2
nb = p+1
nodal_pts, _ = LegendreGaussLobatto(nb)
u = np.ones_like(nodal_pts)

# grad vandermonde
V  = Vandermonde1D(N=nb, x=nodal_pts)
Vr = GradVandermonde1D(N=nb, x=nodal_pts)
Vi = np.linalg.inv(V)
Dr = np.dot(Vr, Vi)

duh = np.dot(Dr, u)
print(np.max(np.abs(duh)))
```

    1.1102230246251565e-16


#### exponential convergence in $p$ for smooth functions

If the function of which we are attempting to take the discrete derivative is smooth (in the sense of infinitely differentiable), then we will see exponential convergence to the analytical solution w/r/t polynomial order of the nodal basis.


```python
ps = [1, 2, 3, 4, 8, 12, 16, 18, 24, 32]
errs = []
fig, ax = plt.subplots(1, 3, figsize=(8, 3))
for p in ps:
    nb = p+1
    nodal_pts, _ = LegendreGaussLobatto(nb)
    view_pts = np.linspace(-1, 1, 100)

    # grad vandermonde
    V  = Vandermonde1D(N=nb, x=nodal_pts)
    Vr = GradVandermonde1D(N=nb, x=nodal_pts)
    Vi = np.linalg.inv(V)
    Dr = np.dot(Vr, Vi)

    # nodal shap
    V = Vandermonde1D(N=nb, x=nodal_pts)
    Vti = np.linalg.inv(V.T)
    P = P_tilde(r=view_pts, N=nb)
    view_shap = np.dot(Vti, P)
    
    u = np.sin(nodal_pts-np.pi/4.)
    du = np.cos(view_pts-np.pi/4)
    duh = np.dot(Dr, u)

    view_duh = np.dot(duh, view_shap)
    err = np.max(np.abs(view_duh - du))
    errs.append(err)
    
    # plot first few
    if p < 4:
        ax[p-1].plot(view_pts, np.sin(view_pts), label=r'$u$')
        ax[p-1].plot(view_pts, du, label=r'$u^\prime$')
        ax[p-1].plot(view_pts, view_duh, '--', label=r'$du_h$')  
        ax[p-1].set_title(r'$p={}$'.format(p))
        ax[p-1].legend()
plt.show()
```


    

    



```python
fig, ax = plt.subplots()
ax.semilogy(ps, errs)
ax.set_xticks(ps)
ax.set_ylabel(r'$||du - du_h||_{L_\infty}$')
ax.set_xlabel('polynomial order p')
ax.set_title('exponential convergence of discrete derivative')
plt.show()
```


    

    


# creation of a 1D master element

Define a simple nodal basis object for 1D problems -- no need to be fancy, pretty much the only thing we need this for is to get nodal shape functions and their derivatives efficiently. The underlying orthonormal Legendre polynomial basis is hidden to the user.


```python
class NodalBasis1D(object):
    """ minimalist nodal basis object: 
    efficiently computes shape functions and their derivatives 
    """
    def __init__(self, p, node_spacing='GAUSS_LOBATTO'):
        self.nb = p + 1
        if node_spacing == 'GAUSS_LOBATTO':
            self.nodal_pts, _ = LegendreGaussLobatto(self.nb)
        elif node_spacing == 'EQUIDISTANT':
            self.nodal_pts = np.linspace(-1, 1, self.nb)
        else: raise ValueError('node_spacing {} not recognized'.format(node_spacing))
    
    def shape_functions_at_pts(self, pts):
        """ computes shape functions evaluated at pts on [-1, 1]
        @retval shap (len(pts), nb) phi_j(pts[i])
        """
        V = Vandermonde1D(N=self.nb, x=self.nodal_pts)
        VTinv = np.linalg.inv(V.T)
        P = P_tilde(r=pts, N=self.nb)
        shap = np.dot(VTinv, P)
        return shap.T
    
    def shape_function_derivatives_at_pts(self, pts):
        """ computes shape function derivatives w/r/t x on [-1, 1]
        @retval shap_der, (Dr in Hesthaven), (len(pts), nb) d/dx phi_j(pts[i])
        """
        V  = Vandermonde1D(N=self.nb, x=self.nodal_pts)
        Vx = GradVandermonde1D(N=self.nb, x=pts)
        Vinv = np.linalg.inv(V)
        shap_der = np.dot(Vx, Vinv)
        return shap_der
```

Define a 1D master element, which is built on top of the 1D basis. 
- Precompute shape functions at the nodal points and Gauss Legendre quadrature points, both are useful for different types of schemes. We use Gauss Legendre points instead of Gauss Lobatto points because they can integrate degree $2n-1$ polynomials exactly instead of $2n - 3$, where $n$ is the number of integration points. We would like to integrate $(\phi_i, \phi_j)_{\hat{K}}$, which is order 2$p$, so to integrate the mass matrix exactly, we need $2p + 1$ points, and common practice is $2p+2$. Since quadrature in 1D is cheap, we opt for the latter.
- Precompute mass matrix $M_{ij} = (\phi_i, \phi_j)$ and stiffness matrices $S_{ij} = \left(\phi_i, \frac{d\phi_j}{dx}\right)$, $K_{ij} = \left(\frac{d\phi_i}{dx}, \frac{d\phi_j}{dx}\right)$. Additionally, store $M^{-1}$, as it is commonly used. Although Hesthaven's method for mass and stiffness matrices are elegant, they rely on the underlying choice of an orthanormal modal basis. Since this class could be overloaded to work with other choices of basis, better to simply compute these matrices with quadrature.


```python
class Master1D(object):
    """ minimalist 1D master object """
    def __init__(self, p, nquad_pts=None, *args, **kwargs):
        self.p, self.nb = p, p+1
        self.basis = NodalBasis1D(p=p, **kwargs)
        self.nodal_pts = self.basis.nodal_pts
        self.nq = 2*self.p + 2 if nquad_pts is None else nquad_pts
        self.quad_pts, self.wghts = GaussLegendre(self.nq)
        
        # shape functions at nodal and quadrature points
        self.shap_nodal, self.dshap_nodal = self.mk_shap_and_dshap_at_pts(self.nodal_pts)
        self.shap_quad,  self.dshap_quad = self.mk_shap_and_dshap_at_pts(self.quad_pts)
        
        # mass, stiffness matrices
        self.M, self.S, self.K = self.mk_M(), self.mk_S(), self.mk_K()
        self.Minv = np.linalg.inv(self.M)
        
        # lifting permuatation matrix L (0s, 1s)
        self.L = self.mk_L()
        
    def mk_shap_and_dshap_at_pts(self, pts):
        shap = self.basis.shape_functions_at_pts(pts)
        dshap = self.basis.shape_function_derivatives_at_pts(pts)
        return shap, dshap
        
    def mk_M(self):
        shapw = np.dot(np.diag(self.wghts), self.shap_quad)
        M = np.dot(self.shap_quad.T, shapw)
        return M
    
    def mk_S(self):
        dshapw = np.dot(np.diag(self.wghts), self.dshap_quad)
        S = np.dot(self.shap_quad.T, dshapw)
        return S
    
    def mk_K(self):
        dshapw = np.dot(np.diag(self.wghts), self.dshap_quad)
        K = np.dot(self.dshap_quad.T, dshapw)
        return K
    
    def mk_L(self):
        L = np.zeros((self.nb, 2))
        L[0, 0]  = 1
        L[-1, 1] = 1
        return L 
    
    @property
    def shap_der(self): 
        """ return the shape derivatives for apps expecting 2, 3D"""
        return [self.dshap_quad]
```

# creation of 1D mesh and DOF handler

## 1D mappings

For 1D problems, the mapping from the master element to physical space elements is somewhat trivial, since there's no reason for the transformation to be anything except affine. Note though, that when the 1D elements are embedded in 2D, then the transformations may be non-affine, in which case we must handle isoparametric mappings and the like. We defer this until later. For an affine mapping, we have the simple mapping

$$x(\xi) = x_L^k + \frac{1 + \xi}{2}(x_R^k - x_L^k)$$

With which we can move the nodal master points to their physical space coordinates.


```python
# build T and P arrays
P = np.linspace(2, 4, 5)
```


```python
class Mesh1D(object):
    def __init__(self, P):
        """ @param P  vertex points, sorted by x position """
        self.verts = P
        self.nElm, self.nEdges = len(self.verts) - 1, len(self.verts)
        self.connectivity = self.build_T()
        connected_one_side = np.bincount(self.connectivity.ravel()) == 1
        self.boundary_verts = np.where(connected_one_side)[0]
        
    def build_T(self):
        """ element connectivity array from 1D vertex list """
        T = np.zeros((self.nElm, 2), dtype=int)
        T[:,0] = np.arange(self.nElm)
        T[:,1] = np.arange(self.nElm) + 1
        return T
```


```python
class dofh_1D(object): pass

class DG_dofh_1D(dofh_1D):
    def __init__(self, mesh, master):
        self.mesh, self.master = mesh, master
        self.n_dof = self.master.nb * self.mesh.nElm
        self.dgnodes = self.mk_dgnodes()
        self.lg = self.mk_lg()
        self.lg_PM = self.mk_minus_plus_lg()
        self.nb, self.nElm = self.master.nb, self.mesh.nElm
        self.ed2elm = self.mk_ed2elm()
    
    def mk_dgnodes(self):
        """ map master nodal pts to element vertices def'd in self.mesh """
        dgn = np.zeros((self.master.nb, self.mesh.nElm))
        master_nodal_pts = np.squeeze(self.master.nodal_pts)
        for elm, elm_verts in enumerate(self.mesh.connectivity):
            elm_vert_pts = self.mesh.verts[elm_verts]
            elm_width = elm_vert_pts[1] - elm_vert_pts[0]
            mapped_pts = elm_vert_pts[0] + (1+master_nodal_pts)/2.*(elm_width)
            dgn[:, elm] = mapped_pts 
        return dgn
    
    def mk_lg(self):
        """ number all dof sequentially by dgnodes """
        node_numbers = np.arange(np.size(self.dgnodes))
        lg = node_numbers.reshape(self.dgnodes.shape, order='F')
        return lg
    
    def mk_minus_plus_lg(self):
        """ (-) denotes element interior, (+) denotes exterior"""
        lg_PM = dict()
        lg_PM['-'] = self.lg[[0, -1], :].ravel(order='F')
        lgP = self.lg[[0, -1],:]
        lgP[0, 1: ] -= 1 # shift nodes to left of first
        lgP[1, :-1] += 1 # shift nodes to right of last
        lg_PM['+'] = lgP.ravel(order='F')
        return lg_PM     
    
    def mk_ed2elm(self):
        """ internal map holding the indicies to reshape vector of values on faces to 
        element edge space (2, nElm), duplicating the values on either side of interior faces
        """
        f2elm = np.zeros((2, self.nElm))
        faces = np.arange(self.mesh.nEdges)
        # numpy magic is doing the following:
        #
        # [[0, 1, 2, 3]
        #  [0, 1, 2, 3]]  - ravel('F') -> [0, 0, 1, 1, 2, 2, 3, 3]
        # 
        # close, but ends duplicated. => trim the ends and reshape to f2elm shape
        # 
        # [[0, 1, 2]
        #  [1, 2, 3]]
        #
        f2elm = np.vstack((faces, faces)).ravel(
            order='F')[1:-1].reshape(f2elm.shape, order='F')
        return f2elm
    
    def edge2elm_ed(self, arr):
        """ internal method to move edge values (defined on the interfaces)
        to values on the "element edge space", the edge dof interior to each element
        @param arr  array formatted on edge space (nFaces,)
        @retval elmEdArr  array formatted on "element edge space" (2, nElm)
        """
        return arr[self.ed2elm]
```

# computation of fluxes

The 'back end' of an explicit DG computation is the unrolled vector of all the problem unknowns. The front end that we'd like to interact with is the dgnodes data structure


```python
def plot_solution(ax, u, dofh):
    """ u formatted like dgnodes """
    for elm in range(dofh.nElm):
        nodal_pts = dofh.dgnodes[:, elm]
        nodal_values = u[:, elm]
        ax.plot(nodal_pts, nodal_values)
    return ax
```


```python
# Low storage Runge-Kutta coefficients LSERK
rk4a = np.array([
         0.0,
        -567301805773.0/1357537059087.0,
        -2404267990393.0/2016746695238.0,
        -3550918686646.0/2091501179385.0,
        -1275806237668.0/842570457699.0])

rk4b = np.array([ 
         1432997174477.0/9575080441755.0,
         5161836677717.0/13612068292357.0,
         1720146321549.0/2090206949498.0,
         3134564353537.0/4481467310338.0,
         2277821191437.0/14882151754819.0])

rk4c = np.array([
         0.0,
         1432997174477.0/9575080441755.0,
         2526269341429.0/6820363962896.0,
         2006345519317.0/3224310063776.0,
         2802321613138.0/2924317926251.0])
```


```python
# constants
π = np.pi

# geometry set up
P = np.linspace(0, 2*π, 10)
mesh1d = Mesh1D(P)
master = Master1D(p=2)
dofh = DG_dofh_1D(mesh1d, master)
mapdgn = np.zeros((dofh.dgnodes.shape[0], 1, dofh.dgnodes.shape[1]))
mapdgn[:,0,:] = dofh.dgnodes
_map = fem_map.Affine_Mapping(master=[master], dgnodes=[mapdgn])
```

We choose numerical fluxes of the form

$$\widehat{au} =\left\{\!\!\left\{au\right\}\!\!\right\} + (1-\alpha)\frac{|a|}{2}  \left[\!\!\left[u\right]\!\!\right]$$

Where $\alpha = 0$ represents an upwinded flux and $\alpha=1$ represents a central flux. These are shown in Hesthaven to be stable for the equation we are interested in solving.


```python
def compute_interior_flux(u, norm, dofh, α):
    """ computes the numerical flux at all of the element interfaces
    @param u  the current solution u, unrolled to a vector
    NOTE: boundary interfaces will be filled with garbage, and must be corrected
    """
    pm = dofh.lg_PM
    
    # equivalent to the flux
    #  \hat{au} = {{au}} + (1-α) * |a|/2 * [[u]]
    # at element interfaces. First and last interface will have garbage.
    flux = a/2*(u[pm['-']] + u[pm['+']]) + (1-α)*np.abs(norm*a)/2.*(u[pm['+']] - u[pm['-']])
    return flux
```

# semi-discrete scheme

Considering the "weak" DG-FEM form,  we have the semi-discrete element local equation

\begin{align}
    \int_K \frac{\partial u_h}{\partial t} v \, dK
    +\int_{K} (au_h) \frac{\partial v}{\partial x} \, dK = 
    -\int_{\partial K} \hat{n}\cdot \widehat{au} v \, d\partial K
\end{align}

Choosing a representation $u=u_i\phi_i$ piecewise polynomial over each element, and the same test space, we have, for a given choice of numerical flux $\widehat{au}$, and noting that in 1D, the normal vectors are simply (-1, +1):

\begin{align}
    \int_K \frac{\partial}{\partial t} (u_i(t) \phi_i(x)) \phi_j(x) \, dx
    +\int_{K} a(u_i(t)\phi_i(x)) \frac{\partial \phi_j(x)}{\partial x} \, dx = 
    -(\widehat{au}(x_R) - \widehat{au}(x_L))
\end{align}

transforming the integrals to the reference element:

\begin{align}
    \int_{\hat{K}} \frac{\partial}{\partial t} (u_i(t) \phi_i(\xi)) \phi_j(\xi) \,|det(J)|\, d\xi
    +\int_{\hat{K}} a(u_i(t)\phi_i(\xi)) \frac{\partial \phi_j(\xi)}{\partial \xi} \, |det(J)|\, d\xi = 
    -(\widehat{au}(x_R) - \widehat{au}(x_L))
\end{align}

This completes the description of the semi-discrete scheme, and we have a choice as to how to compute these integrals. The important part is that since the coefficients $u_i$ vary in time but are constants with respect to space, we can write 

\begin{align}
    &\frac{\partial u_i(t)}{\partial t} \int_{\hat{K}} \phi_i(\xi) \phi_j(\xi) \,|det(J)|\, d\xi
    +au_i\int_{\hat{K}} \phi_i(\xi) \left(\frac{d\xi}{dx}\right)\frac{\partial \phi_j(\xi)}{\partial \xi} \, |det(J)|\, d\xi = 
    -(\widehat{au}(x_R) - \widehat{au}(x_L)) \\
   &\Rightarrow M_K \vec{\frac{du_h}{dt}} + a  S_K \vec{u_h} 
   = - L\, (\widehat{au}(x_R) - \widehat{au}(x_L))
\end{align}

Where we have computed $M_K$ and $S_K$, the mass and stiffness matrices for element $K$. Although we would normally do this with a quadrature rule, we can take advantage of the fact that in 1D (and indeed under any affine mapping from reference to physical element), $J^{-1}$ and $|\det(J)|$ will be constant over the entire element (also note that in 1D, $J^{-1}$ is a 1x1 matrix)<sup>1</sup>. In that case, we can treat both as constants, precompute $M_{\hat{K}}, S_{\hat{K}}$, and multiply the entire element-local equation by $M^{-1}$, giving

\begin{align}
\vec{\frac{du_h}{dt}} &= - a \frac{\det(J)_K}{\det(J)_K}\, J^{-1}_K M_{\hat{K}}^{-1}S^T_{\hat{K}} \vec{u_h} 
  - \frac{1}{\det(J)_K} M^{-1}_K L\, (\widehat{au}(x_R) - \widehat{au}(x_L)) \\
  &= - a \, J^{-1}_K M_{\hat{K}}^{-1}S^T_{\hat{K}} \vec{u_h} 
  - \frac{1}{\det(J)_K} M^{-1}_K L\, (\widehat{au}(x_R) - \widehat{au}(x_L))
\end{align}

Which is a good form for a black box integrator, since we have a "naked" $\frac{du_h}{dt}$, and because the scheme is explicit.

note<sup>1</sup>: $J, J^{-1}$ are 1x1 matrices, and $\det{J}$ is simply $J_{11}$; $J^{-1} = 1/J_{11}$. It's important for the clarity of explicit schemes to understand where these cancellations occur.


```python
def advect_rhs_1D(u, t_local, a, dofh, _map, master, flux_fn, gD, norm):
    
    return u
```


```python
# final time
T = 10

# compute time step size, irrelevant for backward euler
CFL = 0.75
Δx = dofh.dgnodes[1,0] - dofh.dgnodes[0,0]
Δt = CFL/(2*π)*Δx
Δt = Δt / 2

# number of timesteps needed
steps = int(np.ceil(T/Δt))

# initial condition, advection speed
solution = np.zeros((steps, *dofh.dgnodes.shape))
a = 2 * np.pi
solution[0,::] = np.sin(dofh.dgnodes)

LSERK_stages = [0, 1, 2, 3, 4]
t = 0

gD = lambda t: -np.sin(a*t)

# normal vectors, all positive
norm = np.ones((2, dofh.nElm))
norm[0,:] *= -1

# function pointer to something that can compute fluxes
flux_fn = compute_interior_flux

# time loop
RK_resid = np.zeros_like(dofh.dgnodes)
for tstep in range(3):
    u = solution[tstep, ::]
    for s in LSERK_stages:
        t_local = t + rk4c[s]*Δt
        rhsu = advect_rhs_1D(u, t_local, a, dofh, _map, master, flux_fn, gD, norm)
        RK_resid = rk4a[s]*RK_resid + Δt*rhsu
        u += rk4b[s]*RK_resid
    t += Δt
```


```python
u0 = solution[0,:,:]
fix, ax = plt.subplots()
ax = plot_solution(ax, u0, dofh)
```


    

    



```python
pm = dofh.lg_PM
u = u0.ravel()

# normal vectors on interfaces, all positive
norm_faces = np.ones(pm['-'].shape[0])

α = 0

# compute interior fluxes
flux = compute_interior_flux(u, norm_faces, dofh, α)

# compute boundary fluxes
flux[0] = gD(t_local)
flux[-1] = flux[0]

dofh.edge2elm_ed(flux)
```




    array([[-0.37763258,  5.55594966,  6.31799283, -5.55594966, -2.19421588,
             2.19421588, -4.12377695,  4.12377695,  6.31799283],
           [ 5.55594966,  6.31799283, -5.55594966, -2.19421588,  2.19421588,
            -4.12377695,  4.12377695,  6.31799283, -6.31799283]])




```python

```

c50fe9ba54de815b73474e3d3262be8283558b06
ipynb
Jupyter Notebook
codes/Hesthaven_1D.ipynb
CorbinFoucart/FEMexperiment
9bad34d9ed7cbdd740e3a4b67f433779dd53b264
MIT
2018-05-26T22:09:32.000Z
2018-06-25T21:46:32.000Z
codes/Hesthaven_1D.ipynb
CorbinFoucart/FEMexperiment
9bad34d9ed7cbdd740e3a4b67f433779dd53b264
MIT
2018-05-17T21:38:44.000Z
2022-03-11T23:21:25.000Z
codes/Hesthaven_1D.ipynb
CorbinFoucart/FEMexperiment
9bad34d9ed7cbdd740e3a4b67f433779dd53b264
MIT
Qwen/Qwen-72B
1. YES
2. YES