Coding isoparametric mapping#
In this notebook an implementation is provided of numerical integration with isoparametric mapping. At the end some tasks are defined for exploring how isoparametric mapping works. Try to look at the code and understands what happens at the various steps. Go back over to the theory provided in the earlier pages if you cannot follow along.
After you have been over the whole page, you can test yourself by building the function numerical_integration with the functions implemented.
In this exercise the following codeblocks are used as a basis for the function numerical_integration:
the physical coordinates of one of the integration points
the shape function derivatives, first \(\frac{dN}{dx}\) and \(\frac{dN}{dy}\)
the jacobian matrix
the stiffness matrix for 2D poisson, using codes from previous blocks
Then an example is shown using the function numerical_integration for a quadrilateral element in 2D for Poisson. This element is also plotted using the physical and reference coordinates.
After you have seen the implementation, it is up to you to answer some questions. The solutions are provided as well.
\(\newcommand{\pder}[2]{\frac{\partial #1}{\partial #2}}\) \(\newcommand{\hpder}[2]{\displaystyle\frac{\partial #1}{\partial #2}}\)
# Import dependencies
import numpy as np
import matplotlib.pyplot as plt
Physical coordinates#
For two natural coordinates \(\xi\) and \(\eta\), the reference quad is defined over the domain \((\xi\in[-1,1],\eta\in[-1,1])\) and its physical counterpart is arbitrarily positioned in 2D space with nodal coordinates \((x_i,y_i),\ i=1\ldots4\). Shape functions are defined for the reference element as:
The \((\xi,\eta)\)-coordinates of any point inside the reference element are mapped to \((x,y)\)-coordinates in physical space through a mapping:
The implementation of this can be seen in the code below. Node that the coordinates are first extracted from a coordinate list using a coords function.
def coords(coordinates):
# Nodal coordinates of the quadrilateral element
x1, y1 = coordinates[0]
x2, y2 = coordinates[1]
x3, y3 = coordinates[2]
x4, y4 = coordinates[3]
x_coords = np.array([x1, x2, x3, x4])
y_coords = np.array([y1, y2, y3, y4])
return x_coords, y_coords
def physical_coordinates(coordinates, xi, eta):
# Extract coordinates
x_coords, y_coords = coords(coordinates)
# Shape functions for a quadrilateral element
N_1 = 1/4 * (1 - xi) * (1 - eta)
N_2 = 1/4 * (1 + xi) * (1 - eta)
N_3 = 1/4 * (1 + xi) * (1 + eta)
N_4 = 1/4 * (1 - xi) * (1 + eta)
# Calculate physical coordinates (x, y)
x = N_1 * x_coords[0] + N_2 * x_coords[1] + N_3 * x_coords[2] + N_4 * x_coords[3]
y = N_1 * y_coords[0] + N_2 * y_coords[1] + N_3 * y_coords[2] + N_4 * y_coords[3]
return x, y
\(\renewcommand{\pder}[2]{\frac{\partial #1}{\partial #2}}\) \(\renewcommand{\hpder}[2]{\displaystyle\frac{\partial #1}{\partial #2}}\)
Derivatives shape functions#
In order to create the Jacobian matrix, we first need to create the derivatives of the shape functions with respect to \(\xi\) and \(\eta\), which are \(\pder{N_i}{\xi}\) and \(\pder{N_i}{\eta} \).
def shape_functions_derivatives(xi, eta):
# Derivatives of shape functions with respect to xi
dN_dxi_1 = -1/4 * (1 - eta)
dN_dxi_2 = 1/4 * (1 - eta)
dN_dxi_3 = 1/4 * (1 + eta)
dN_dxi_4 = -1/4 * (1 + eta)
# Derivatives of shape functions with respect to eta
dN_deta_1 = -1/4 * (1 - xi)
dN_deta_2 = -1/4 * (1 + xi)
dN_deta_3 = 1/4 * (1 + xi)
dN_deta_4 = 1/4 * (1 - xi)
ref_derivatives = np.array([[dN_dxi_1, dN_dxi_2, dN_dxi_3, dN_dxi_4],
[dN_deta_1, dN_deta_2, dN_deta_3, dN_deta_4]])
return ref_derivatives
\(\renewcommand{\pder}[2]{\frac{\partial #1}{\partial #2}}\) \(\renewcommand{\hpder}[2]{\displaystyle\frac{\partial #1}{\partial #2}}\)
Jacobian#
The Jacobian imatrix can now be defined as derivatives of physical coordinates with respect to natural coordinates:
Note that not only \(J\) but also the determinant of \(J\) is returned.
def Jacobian(coordinates, ref_derivatives):
# Extract coordinates
x_coords, y_coords = coords(coordinates)
J = np.array([
[np.dot(ref_derivatives[0,:], x_coords),
np.dot(ref_derivatives[0,:], y_coords)],
[np.dot(ref_derivatives[1,:], x_coords),
np.dot(ref_derivatives[1,:], y_coords)]
])
return J
\(\renewcommand{\pder}[2]{\frac{\partial #1}{\partial #2}}\) \(\renewcommand{\hpder}[2]{\displaystyle\frac{\partial #1}{\partial #2}}\)
We also would like to have the shape function derivatives with respect to \(x\) and \(y\). Thefore we need the inverse of the Jacobian matrix \(J\).
def shape_functions_xy(J, ref_derivatives):
# Calculate the inverse of the Jacobian matrix
J_inv = np.linalg.inv(J)
# Derivatives of shape functions with respect to x and y
dN_dx = J_inv[0, 0] * ref_derivatives[0,:] + \
J_inv[0, 1] * ref_derivatives[1,:]
dN_dy = J_inv[1, 0] * ref_derivatives[0,:] + \
J_inv[1, 1] * ref_derivatives[1,:]
return dN_dx, dN_dy
\(\renewcommand{\pder}[2]{\frac{\partial #1}{\partial #2}}\) \(\renewcommand{\hpder}[2]{\displaystyle\frac{\partial #1}{\partial #2}}\)
Poisson in 2D#
For Poisson it has been learned that
in which:
where:
For the element matrix \(\mathbf{K}_e\) the matrices need to be evaluated at the different integration points and then summed whilst incorporating their different weights:
def B_element_matrix(dN_dx, dN_dy):
B = np.zeros((2, len(dN_dx)))
for i in range(len(dN_dx)):
B[0,i] = dN_dx[i]
B[1,i] = dN_dy[i]
return B
def K_element_matrix(B, nu, J):
detJ = np.linalg.det(J)
K = nu * np.dot(B.T, B) * detJ
return K
Quadrilateral element in 2D for Poisson#
All the steps necessary to get to the element stiffness matrix \( \mathbf{K}_{e} \) matrix can now be used. The only thing left for us to do is to combine it all in a logical order in a function. This has been implemented below in the function numerical_integration. The function needs as input the physical coordinates of the corner points of the element, poisson’s ratio and the number of integration points. It returns the Jacobians at the different integration points, the \(\xi\)- and \(\eta\)-coordinates in the reference configuration of these integration points, and the corresponding weights of the element function. It also returns the physical coordinates of the integration points and the \( \mathbf{K}_{e} \) matrix.
Note
The number specified for \(n_{ip}\) is the number of integration points for each axis. This means that for \(n_{ip} = 2\), we actually have \(2 \cdot 2 = 4\) integration points. And for \(n_{ip} = 7\) we have \(7 \cdot 7 = 49\) integration points.
def numerical_integration(coordinates, nu, n_ip):
x, w = np.polynomial.legendre.leggauss(n_ip)
gauss_pts = np.array(np.meshgrid(x,x,indexing='ij')).reshape(2,-1).T
xi_coords = np.zeros(len(gauss_pts))
eta_coords = np.zeros(len(gauss_pts))
weights = (w*w[:,None]).ravel()
K_element = np.zeros((4,4))
x_coords = np.zeros(len(xi_coords))
y_coords = np.zeros(len(eta_coords))
Jacobians = np.zeros((len(weights), 2, 2))
for i in range(len(gauss_pts)):
xi_coords[i] = gauss_pts[i, 0]
eta_coords[i] = gauss_pts[i, 1]
for i in range(len(weights)):
xi = xi_coords[i]
eta = eta_coords[i]
weight = weights[i]
x_coords[i], y_coords[i] = physical_coordinates(coordinates, xi, eta)
ref_derivatives = shape_functions_derivatives(xi, eta)
J = Jacobian(coordinates, ref_derivatives)
Jacobians[i] = J
dN_dx, dN_dy = shape_functions_xy(J, ref_derivatives)
B_integration_point = B_element_matrix(dN_dx, dN_dy)
K_integration_point = K_element_matrix(B_integration_point, nu, J)
K_element += weight * K_integration_point
return Jacobians, xi_coords, eta_coords, x_coords, y_coords, K_element
We can run the code when we define the input. Consider a distorted quadrilateral element with physical nodal coordinates:
\(\mathbf{x}_1 = (0,0), \quad\mathbf{x}_2=(4,0), \quad\mathbf{x}_3 = (6,2), \quad\mathbf{x}_4 = (0,5)\),
\(\nu = 0.2\), and 2 integration points.
coordinates = [(0, 0), (4, 0), (6, 2), (0, 5)]
nu = 0.2
n_ip = 2
Jacobians, xi_coords, eta_coords, x_coords, y_coords, K_element = numerical_integration(coordinates, nu, n_ip)
Now we can plot the obtained results. Below you can see the results. Note that the boxes around the element are also plotted.
Show code cell source
x_box, y_box = coords(coordinates)
x_box = np.append(x_box, x_box[0])
y_box = np.append(y_box, y_box[0])
xi_box = np.array([-1, -1, 1, 1, -1])
eta_box = np.array([-1, 1, 1, -1, -1])
plt.figure(figsize=(16, 8))
plt.subplot(1,2,1)
plt.plot(x_box, y_box)
plt.plot(x_coords, y_coords, marker='x', linestyle='')
# Annotate each point with a label
for i in range(len(x_coords)):
label = f'{i}: ({x_coords[i]:.2f}, {y_coords[i]:.2f})'
plt.text(x_coords[i], y_coords[i], label, fontsize=8, ha='left')
plt.title('Physical coordinates')
plt.xlabel('x-coordinate')
plt.ylabel('y-coordinate')
plt.subplot(1,2,2)
plt.plot(xi_box, eta_box)
plt.plot(xi_coords, eta_coords, marker='x', linestyle='')
# Annotate each point with a label
for i in range(len(xi_coords)):
label = f'{i}: ({xi_coords[i]:.2f}, {eta_coords[i]:.2f})'
plt.text(xi_coords[i], eta_coords[i], label, fontsize=8, ha='left')
plt.title('Reference coordinates')
plt.xlabel('\u03BE-coordinate')
plt.ylabel('\u03B7-coordinate');
Task 1: Jacobian
The Jacobian matrix below has been found in a previous exercise (see Isoparametric mapping) as a function of \(\xi\) and \(\eta\). Print the Jacobian of the left bottom coordinate (index 0) and fill in the matrix below. Check if it matches.
# YOUR_CODE_HERE
Solution
The following code is an example how to do so:
print(f'The function will give the following Jacobian:')
print(Jacobians[0])
xi = -0.58
eta = -0.58
print(f'The bottom left point ({xi}, {eta}) will give the following Jacobian:')
J = np.array([[1/2*(5+eta), -3/4*(1+eta)],
[1/2*(1+xi), 1/4*(7-xi)]])
print(J)
Task 2: Weights for numerical integration
For integration in 2D on a reference element, the weights at which the numerical integration has been performed should sum up to a certain value. For a 1D reference element, the sum of weights is equal to 2. What number do you expect for 2D? Can you reason why? First use 3 and then 8 integration points (n_ip). Does this matter for the sum of the weights? Print the sum of the weights for both cases and come up with a reasonable explanation. Rewatch the YouTube video on numerical integration if you have to.
# YOUR_CODE_HERE
Solution
For 2D the weights should sum up to 4. This can be explained by the integral below, that should equal the area of the reference element (which is 2x2, as both \(\xi\) and \(\eta\) are in [-1, 1]) Note that the number of integration points does not matter for the individual weights.
The following lines of code will show this:
n_ip = 3
x, w = np.polynomial.legendre.leggauss(n_ip)
weights = (w*w[:,None]).ravel()
print(f'Sum weights = {np.sum(weights):.3f}')
Task 3: Differences Jacobian for rectangular and distorted elements
Take a rectangular elements, for example:
\(\mathbf{x}_1 = (0,0), \quad\mathbf{x}_2=(4,0), \quad\mathbf{x}_3 = (4,4), \quad\mathbf{x}_4 = (0,4)\)
Now take a distorted element, for example:
\(\mathbf{x}_1 = (0,0), \quad\mathbf{x}_2=(6,0), \quad\mathbf{x}_3 = (20,4), \quad\mathbf{x}_4 = (10,4)\)
If you increase the number of integration points, the Jacobian will not change for rectangular elements. However, for distorted elements, the Jacobian will change. Can you explain why this is the case? Use the numerical_integration
function multiple times to output the results for both cases, in which you vary the number of integration points. Hint: think of the mapping between \((\xi, \eta)\) and \((x, y)\) and its implications for the Jacobian, which is a set of derivatives.
#YOUR_CODE_HERE
Solution
You can use the code below in order to look at the results and see that the above statement is true. For rectangular elements, the mapping between the local coordinates (natural coordinates) and global coordinates is linear. This means that the Jacobian matrix, which contains the partial derivatives of the global coordinates with respect to the local coordinates, is constant across the element. For distorted elements, the mapping between the local and global coordinates is generally non-linear. The shape functions used in isoparametric elements describe a more complex relationship, which can vary within the element.
If you increase the number of integration points in an element, the location of the integration points will change. Since for rectangular elements the Jacobian is constant and irrespective of the integration point location across the element, it will not change. Since the mapping is non-linear, the Jacobian changes depending on the location within the element.
print('Rectangular element')
coordinates = [(0, 0), (4, 0), (4, 4), (0, 4)]
nu = 0.2
n_ip_list = np.array([1, 2, 3, 4])
for i in range(len(n_ip_list)):
n_ip = n_ip_list[i]
Jacobians, xi_coords, eta_coords, x_coords, y_coords, K_element = numerical_integration(coordinates, nu, n_ip)
print(n_ip)
print(Jacobians[0])
print('Distorted element')
coordinates_dist = [(0, 0), (6, 0), (20, 4), (10, 4)]
for i in range(len(n_ip_list)):
n_ip = n_ip_list[i]
Jacobians, xi_coords, eta_coords, x_coords, y_coords, K_element = numerical_integration(coordinates_dist, nu, n_ip)
print(n_ip)
print(Jacobians[0])
Task 4: Build the numerical integration function yourself
To get you started, some code has already been written. Good luck!
def numerical_integration(coordinates, nu, n_ip):
x, w = np.polynomial.legendre.leggauss(n_ip)
gauss_pts = np.array(np.meshgrid(x,x,indexing='ij')).reshape(2,-1).T
weights = (w*w[:,None]).ravel()
# --------------
# YOUR_CODE_HERE
# --------------
for i in range(len(gauss_pts)):
# --------------
# YOUR_CODE_HERE
# --------------
for i in range(len(weights)):
# --------------
# YOUR_CODE_HERE
# --------------
return Jacobians, xi_coords, eta_coords, x_coords, y_coords, K_element