# HDG 1D

If you are reading this introduction outside of the notebook, note that python-only and notebook versions of the code are available. Note that within the notebook, these previous links are not available, you are already all set to start this hands-on session. You may want to have a look at this jupyter tutorial if you are not familiar with jupyter.

## 1. Imports

import numpy as np
import random
import numpy.polynomial.polynomial as polynom
import numpy.linalg as linalg
from scipy import sparse
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt


## 2. Mesh Class

class Mesh:
def __init__(self, omega, number_media, length_media, rho_media, mu_media,
h_glob, order_glob):
#omega : frequency
#number_media : number of layers
#length_media : a vector of size number_media containing the length of each layer
#rho_media : a vector of size number_media containing the density of each layer
#mu_media : a vector of size number_media containing the bulk modulus of each layer
length = 0  #the total length
rho = []  #a vector containing the density of each cell of the mesh
mu = []  #a vector containing the bulk modulus of each cell of the mesh
order = [
]  #a vector containing the polynomial order of each cell of the mesh
nbpoints = 1
# the number of nodes of the mesh
nbele = 0
# the number of elements of the mesh
coor = [
0
]  #a vector of size nbpoints, containing the coordinates of each node of the mesh
ele = [
]  # a table of size nbele x 2 : the first columns gives the first node of the celle, the second gives the second nodes. It is useless in 1D, but we set it in order to mimic 2D and 3D
L = 0
for i in range(0, number_media):
#we loop over all the layers
L = L + length_media[
i]  #we increase the total length by the size of the current layer
c = np.sqrt(
mu_media[i] /
rho_media[i])  #we compute the speed of the acoustic wave
h = h_glob  #the size of the cells in the layer #TBC (exo 8)
order_cell = order_glob  #the order of the cells in the layer
N = int(np.floor(length_media[i] /
h))  #we compute the number of cells of the layer
for j in range(0, N - 1):
#we construct the cells of the layer i
nbele = nbele + 1  #we increase the number  of cells
nbpoints = nbpoints + 1  #we increase the number  of nodes
length = length + h  #we increase the length by h
coor.append(length)  #this is the coordinate of the new node
rho.append(
rho_media[i])  #the density of the cell is rho_media[i]
mu.append(
mu_media[i])  #the bulk modulus of the cell is mu_media[i]
order.append(
order_cell)  #the bulk modulus of the cell is order_global
ele.append(
nbpoints - 2
)  #the two nodes of the elements are nbpoints-2 and nbpoints-1
ele.append(nbpoints - 1)
#now, we define the last cell of the layer
nbele = nbele + 1
nbpoints = nbpoints + 1
length = L
coor.append(length)
rho.append(rho_media[i])
mu.append(mu_media[i])
order.append(order_cell)
ele.append(nbpoints - 2)
ele.append(nbpoints - 1)

self.ele = np.reshape(np.array(ele), [nbele, 2])
self.coor = np.array(coor)
self.rho = np.array(rho)
self.mu = np.array(mu)
self.nbele = nbele
self.nbpoints = nbpoints
self.order = order


## 3. Basis Class

class Basis:
def __init__(self, p):
#p is the order of the basis
#we build the p+1 elementary polynomials of degree 1
poly_elem = np.zeros(
(p + 1, 2)
)  #poly_elem is table (p+1)x2 : each line contains the coefficients a and b such that the polynom reads as a+bx.
for i in range(0, p + 1):
#the root of polynonimal i is i/p, so that a=-i/p and b=1;
poly_elem[i, 0] = -(i) / p
poly_elem[i, 1] = 1
#we initialize the basis
phi = []
#phi is a vector of vector : phi[i] is the basis function i which is a polynomial of degree p, it represented by a vector of size p+1 containing the coefficients of the polynomial fom the lowest to the highest
for i in range(0, p + 1):
#the basis function i is the product of all elementary polynomials except polynomial i. It is normalized so that phi[i](i)=1
phi.append()
for j in range(0, p + 1):
if i != j:
phi[i] = 000000  #TBC (ex 1)
#we multiply phi[i] by all elementary polynomials except polynomial i
phi[i] = phi[i] / 000000  #TBC (ex 1)
#we normalize phi[i]
self.phi = phi
massmat = np.zeros((p + 1, p + 1))  #we initialize the mass matrix
stiffmat = np.zeros(
(p + 1, p + 1))  #we initialize the stiffness matrix
for i in range(0, p + 1):
for j in range(0, p + 1):
massmat[i, j] = 000000  #TBC (ex 2)
stiffmat[i, j] = 000000  #TBC (ex 2)
self.massmat = massmat
self.stiffmat = stiffmat


### 3.1. Action 1: Construction of the basis functions (Class Basis)

#### 3.1.1. Multiplying two polynomials

The python command for multiplying two polynomials q and r is polynom.polymul(q,r). Compute the multiplication of polynomials phi[i] and poly_elem[j,:] and store it in phi[i].

#### 3.1.2. Evaluating a polynomial

The python command for evaluating a polynomial q at point x is polynom.polyval(x,q). Normalize phi[i] by its value at point i/p.

### 3.2. Action 2: Construction of the reference matrices (Class Basis)

#### 3.2.1. Primitive of a polynomial

The python command for computing the primitive of a polynomial P is polynom.polyint(P). Hence, $$\int_0^{1}P(x)dx$$ can be computed thanks to sum(polynom.polyint(P)). Compute: $massmat(i,j) = \int_0^1 \phi_i(x) \phi_j(x) dx$.

#### 3.2.2. Derivative of a polynomial

The python command for computing the derivative of a polynomial is polynom.polyder(P). Compute: $stiffmat(i,j)= \int_0^1 \phi_i(x) \frac{\partial\phi_j(x)}{\partial x} dx$

## 4. Elem Class

#This class contains the local matrix and rhs
class Elem:
def __init__(self, x0, x1, rho, mu, order, alpha, basis, omega, xs):
#x0 : left coordinate of the cell
#x1 :  right coordinate of the cell
#rho : density of the cell
#mu : bulk modulus of the cell
#order : order of the cell
#alpha : penalization coefficient
#basis : the basis functions
#omega : the frequency
#xs : the position of the source, if any
n = order + 1
#the number of degree of freedom in the cell
h = x1 - x0  #the length of the cell
mat = np.zeros(
(2 * n, 2 * n), dtype=complex)  #initialization of the matrix
rhs = np.zeros((2 * n), dtype=complex)  #initialization of the rhs

#the four blocks of the matrix
mat[:n, :n] = 000000  #TBC(exo3)
mat[n:, :n] = 000000  #TBC(exo3)
mat[:n, n:] = 000000  #TBC(exo3)
mat[n:, n:] = 000000  #TBC(exo3)

mat[n, n] = 000000  #TBC(exo3)
mat[2 * n - 1, 2 * n - 1] = 000000  #TBC(exo3)

F = np.zeros((2 * n, 2))  #the right hand side
F[0, 0] = 000000  #TBC(exo3)
F[n - 1, 1] = 000000  #TBC(exo3)
F[n, 0] = 000000  #TBC(exo3)
F[2 * n - 1, 1] = 000000  #TBC(exo3)

#computation of the source if it is inside the element
if (xs >= x0 and xs <= x1):
#we compute the coordinates of the source in the reference element
xs_ref = (xs - x0) / h
for i in range(n):
#we compute the values of basis functions at point xs_ref
rhs[i] = polynom.polyval(xs_ref, basis.phi[i])
#the source is only on the pressure
rhs[n:] = 0
rhs = linalg.solve(mat, rhs)
else:
rhs[:] = 0

self.mat = linalg.solve(mat, F)
self.rhs = rhs


### 4.1. Action 3: Construction of the local matrices and rhs on each element (Class Elem)

#### 4.1.1. Compute the four blocks of the local matrix (without the penalization coefficients)

We have:

1. $$M(0:n-1, 0:n-1) = \frac{i\omega}{\mu} * Mass * h$$
2. $$M(n:2*n-1, 0:n-1) = -Stiff$$
3. $$M(0:n-1, n:2*n-1) = Stiff^T$$
4. $$M(n:2*n-1, n:2*n-1) = i\omega * \rho * Mass * h$$

Note that:

• the $$i$$ imaginary number is defined by 1j in python;
• the mass matrix Mass is in basis.massmat;
• the stiffness matrix Stiff is in basis.stiffmat and its transpose in basis.stiffmat.T.

#### 4.1.2. Add the penalization coefficients

We have:

1. $$M(n,n) = M(n,n)-\alpha$$
2. $$M(2 * n-1, 2*n-1) = M(2*n-1,2*n-1)-\alpha$$

## 5. Class Sparse

class Sparse:
def __init__(self, alpha, elem, nbpoints, order):
#alpha : penalization coefficient
#elem : the local matrices associated to each elements
#nbpoints : the number of vertices of the mesh
#order : a vector containing the order of each element

#the global matrices is stored in coordinate format :
A = []
#the vector containing the non-zero values
I = []
# the vector containing the lines of the non-zero values
J = []
# the vector containing the columns of the non-zero values
#we have Mat(I[i], J[i])=A[i]
B = []
# the Right Hand Side

#Left point of the mesh
p = order

res = 000000  #TBC (exo4)
A.append(res)
I.append(0)
J.append(0)

res = 000000  #TBC (exo4)
A.append(res)
I.append(1)
J.append(0)

#right hand side
res = 000000  #TBC (exo4)
B.append(res)

#now we run on all the vertices of the mesh (except the left and right ones)
for i in range(1, nbpoints - 1):
#we get the order of elem i-1 and i
p1 = order[i - 1]
p2 = order[i]

res = 000000  #TBC (exo4)
A.append(res)
I.append(i - 1)
J.append(i)

res = 000000  #TBC (exo4)
A.append(res)
I.append(i)
J.append(i)

res = 000000  #TBC (exo4)
A.append(res)
I.append(i + 1)
J.append(i)

#right hand side
res = 000000  #TBC (exo4)
B.append(res)

#right point of the mesh
p = order[nbpoints - 2]
res = 000000  #TBC (exo4)
A.append(res)
I.append(nbpoints - 2)
J.append(nbpoints - 1)

res = 000000  #TBC (exo4)
A.append(res)
I.append(nbpoints - 1)
J.append(nbpoints - 1)

res = 000000  #TBC (exo4)
B.append(res)

# Left and right boundary condition. Comment only if there is a point soucre
B = 0
B[nbpoints - 1] = -np.sin(omega * mesh.coor[nbpoints - 1])
self.matrix = sparse.coo_matrix(
(A, (I, J)), shape=(nbpoints, nbpoints))
self.rhs = B


### 5.1. Action 4: Construction of the global matrix (Class Elem) and rhs

#### 5.1.1. For internal vertices

The element left to vertex $$i$$ is $$i-1$$, the element right to vertex $$i$$ is $$i$$. We denote by $$p_1$$ the order of cell $$i-1$$ and by $$p_2$$ the order of cell $$i$$. Check that:

1. $$A(i-1,i)=-M_{i-1}(p_1, 0)-\alpha M_{i-1}(2p_1+1,0)$$
2. $$A(i,i)=-M_{i-1}(p_1,1)-\alpha M_{i-1}(2p_1+1,1)+M_i(0,0)-\alpha M_i(p_2+1,0)+2\alpha$$
3. $$A(i+1,i)=M_i(0,1)-\alpha M_i(p_2+1,1)$$
4. $$B(i)=F_{i-1}(0)-F_i(p_2)+\alpha(F_i(2p_2+1)+F_{i-1}(p1+1))$$

Note that $$M_{i}(j, k)$$ corresponds to elem[i].mat[j,k] and that Fi(j) corresponds to elem[i].rhs[j] in the code.

#### 5.1.2. For left vertex

The element right to vertex $$0$$ is $$0$$, we denote by $$p$$ its order. Check that:

1. $$A(0,0)=M_0(0,0)-2 \alpha M_0(p+1,0)+2\alpha$$
2. $$A(1,0)=M_0(0,1)-2 \alpha M_0(p+1,1)$$
3. $$B(0)=-F_0(p)+2\alpha F_0(2p+1)$$

#### 5.1.3. For right vertex

The element left to vertex $$nbpoints-1$$ is $$nbpoints-2$$, we denote by $$p$$ its order. Check that:

1. $$A(nbpoints-2,nbpoints-1)=-M_{nbpoints-2}(p,0)-2 \alpha M_{nbpoints-2}(2p+1,0)$$
2. $$A(nbpoints-1,nbpoints-1)=-M_{nbpoints-2}(p,1)-2 \alpha M_{nbpoints-2}(2p+1,1)+2\alpha$$
3. $$B((nbpoints-1)=F_{nbpoints-2}(0)+2\alpha F_{nbpoints-2}(p+1)$$

## 6. Main Program

alpha = 0.5  #penalization parameter
omega = 10  #frequency
x_source = 15  #position of the source (if no source, put it outside of the domain)
nb_media = 1  #number of layer
points_plot = 10  #number of points on each cell for plotting the solution
h_glob = 0.1  #Global size of the cells
order_glob = 1  #Global order of the cells
length_media =   #Length of the layers
rho_media =   #density of the layers
mu_media =   #bulk modulus of the layers
mesh = Mesh(omega, nb_media, length_media, rho_media, mu_media, h_glob,
order_glob)  #we create the mesh
basis = list(range(max(mesh.order) + 1))
for p in set(mesh.order):
print("building basis functions of order", p)
basis[p] = Basis(p)  #we define all the basis functions we need

elem = list(range(mesh.nbele + 1))
for i in range(mesh.nbele):
#we run over all the cells of the mesh
node = mesh.ele[i]
x0 = mesh.coor[node]
x1 = mesh.coor[node]
rho = mesh.rho[i]
mu = mesh.mu[i]
order = mesh.order[i]
elem[i] = Elem(x0, x1, rho, mu, order, alpha, basis[order], omega,
x_source)  #we construct the elementary basis functions

linsys = Sparse(alpha, elem, mesh.nbpoints,
mesh.order)  # we create the global linear system
Lambda = 000000  #TBC (exo 5) #we solve the linear system
#plt.plot(np.imag(Lambda))
#plt.show()
#Computation of the solution

#Now, we reconstruct the solution on each element of the mesh
P = []
U = []
X = []
P_ex = []
U_ex = []
err_P = 0  #for computing the error between the approximate solution and the exact one
err_U = 0
L2_P = 0  #L2 norm of the exact solution
L2_U = 0
#we run over all the elements of the mesh
for i in range(mesh.nbele):
#we find the two vertices of the element
node = mesh.ele[i, :]
#we get the two lambda associated to these vertices
lambda_loc = Lambda[node]
#we can now compute locally the solution on the element
W_K = 000000  #TBC (exo 6)
p = mesh.order[i]
#the first degrees of freedom are associated to the velocity
U_K = W_K[p + 1:]
#the last degrees of freedom are associated to the pressure
P_K = W_K[:p + 1]
#we compute the value of P and U on the points of the submesh
for j in range(points_plot + 1):
x = j / points_plot
xcoor = mesh.coor[node] + x * (mesh.coor[node] -
mesh.coor[node])
X.append(xcoor)
P_exact = np.sin(omega * xcoor)
U_exact = -1j * np.cos(omega * xcoor)
P_ex.append(P_exact)
U_ex.append(U_exact)
valP = 0
valU = 0
#We get the order of the cell
for k in range(p + 1):
val = 000000  #TBC (exo 6)
valP = valP + val * P_K[k]
valU = valU + val * U_K[k]
P.append(valP)
U.append(valU)
err_P = err_P + abs(valP - P_exact) * abs(valP - P_exact)
err_U = err_U + abs(valU - U_exact) * abs(valU - U_exact)
L2_P = L2_P + abs(P_exact) * abs(P_exact)
L2_U = L2_U + abs(U_exact) * abs(U_exact)
#computation of relative  L2 error
print("relative L2 error on P :", np.sqrt(err_P / L2_P))
print("relative L2 error on U :", np.sqrt(err_U / L2_U))

fig, axs = plt.subplots(2, 2)
ax1 = axs[0, 0]
ax1.plot(X, np.real(P), 'g')
ax1.plot(X, np.real(P_ex), 'b--')
ax1.set_title('P_real')
ax2 = axs[1, 0]
ax2.plot(X, np.imag(P), 'g')
ax2.plot(X, np.imag(P_ex), 'b--')
ax2.set_title('P_imag')
ax3 = axs[0, 1]
ax3.plot(X, np.real(U), 'g')
ax3.plot(X, np.real(U_ex), 'b--')
ax3.set_title('U_real')
ax4 = axs[1, 1]
ax4.plot(X, np.imag(U), 'g')
ax4.plot(X, np.imag(U_ex), 'b--')
ax4.set_title('U_imag')
plt.show()
fig, axs = plt.subplots(2, 2)
ax1 = axs[0, 0]
ax1.plot(X, np.real(P) - np.real(P_ex), 'g')
ax1.set_title('P_real')
ax2 = axs[1, 0]
ax2.plot(X, np.imag(P) - np.imag(P_ex), 'g')
ax2.set_title('P_imag')
ax3 = axs[0, 1]
ax3.plot(X, np.real(U) - np.real(U_ex), 'g')
ax3.set_title('U_real')
ax4 = axs[1, 1]
ax4.plot(X, np.imag(U) - np.imag(U_ex), 'g')
ax4.set_title('U_imag')
plt.show()


### 6.1. Action 5: Solution to the global system

The solution to a sparse linear system $$Ax=b$$ is given by x=spsolve(A,b). Compute Λ, solution to linsys.matrix \Lambda= linsys.rhs.

### 6.2. Action 6: Reconstruction of the solution, element by element

#### 6.2.1. The matrix vector product is given by np.dot(A,x)

Compute $$W_K=A_K^{-1} \lambda_{loc}+F_K$$ for each element $$K$$ (we remind that AK-1=elem[i].mat and FK is elem[i].rhs)

#### 6.2.2. In order to compute the solution at point x, compute φk(x) for each k

We remind that the value of a polynom p at point x is given by polynom.polyval(x,p)

### 6.3. Action 7: Validation and plotting

#### 6.3.2. Set the order to 1

Check that the relative error is divided by 4 when you divide the space step h by 2

#### 6.3.3. Set the order to p=2…5

Check that the relative error is divided by 2p+1 when you divide the space step h by 2

### 7.1. Action 8: Role of h and p

#### 7.1.1. Set h to 0.1

How does the error behaves when you increase the order, for p from 1 to 10 ?

#### 7.1.2. Choosing the number of point per wavelength at a given order

The wavelength is $$\frac{2\pi c}{\omega}$$, with $$c=\sqrt{\frac{\mu}{\rho}}$$. Determine (roughly) the minimum number of point pers wavelength to obtain a relative error of 0.001, for p from 1 to 5.

#### 7.1.3. Choosing the number of point per wavelength at a given h

Conversely, determine the minimal order to obtain a relative error of 0.001, for a number of point per wavelength from 15 to 1.

Modify the class Mesh in order to define h on each layer, depending on c, omega and orderglob How does the error behaves if you increase the frequency ?

Modify the class Mesh in order to define p on each layer, depending on c, omega and hglob How does the error behaves if you increase the frequency ?

### 7.2. Action 9: Point source

#### 7.2.1. Comment the two lines

• B=0;
• B[nbpoints-1]=-np.sin(omega*mesh.coor[nbpoints-1]);

They are at the end of the Sparse class.

### 7.3. Action 10: Multilayers

Set n_media=3, set length_media=[1,1,1 ], rho_media=[1,2,1] and mu_media=[1,1,2]. Run the code, using either the h- or the p-adaptivity.

## 8. Correction

The corrections below are to be fetched outside from the notebook. If you are in the notebook, you can directly load them from /home/hpcs-readonly/notebook (you can still retrieve the below links at the bottom of this page).

### 8.2. notebook version

Created: 2021-06-10 Thu 18:01

Validate