<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Analytic Finite Element Method
This is a FEM implementation in a CAS and is intended to be useful
for understanding the implementation of the method.
The code solves the Navier-Cauchy equations for elasticity in 2D.
Step1:  Interpolation functions and matrices
Step2:  The interpolation matrix is a matrix with the interpolation
functions reordered for each degree of freedom.
Step3:  The local derivatives matrix is formed with the derivatives of the interpolation functions
Step4:  Gauss integration
Step5:  Local matrices generation
The Jacobian of the transformation is given by the product of the derivative
matrix and the coordinates of each element
Step6:  We can re-arrange the derivatives of the interpolation function as a matrix that
translates from displacements to strains, commonly denoted by $B$
Step7:  With these elements we can form the local stiffness and mass matrices.
Step8:  Assembly process
Step9:  Example
Step10:  Example
Step11:  Example
<ASSISTANT_TASK:>
Python Code:
from sympy import *
init_session()
def const_model(E, nu, const="plane_stress"):
    if const == "plane_stress":
        fac = E/(1 - nu**2)
        C = fac*Matrix([
            [1, nu, 0],
            [nu, 1, 0],
            [0, 0, (1 - nu)/2]])
    elif const == "plane_strain":
        fac = E*(1 - nu)/((1 - 2*nu)*(1 + nu))
        C = fac*Matrix([
            [1, nu/(1 - nu), 0],
            [nu/(1 - nu), 1, 0],
            [0, 0, (1 - 2*nu)/(2*(1 - nu))]])
    return C
r, s = symbols("r s")
Explanation: Analytic Finite Element Method
This is a FEM implementation in a CAS and is intended to be useful
for understanding the implementation of the method.
The code solves the Navier-Cauchy equations for elasticity in 2D.
End of explanation
N = S(1)/4 *Matrix([
[(1 - r)*(1 - s)],
[(1 + r)*(1 - s)],
[(1 + r)*(1 + s)],
[(1 - r)*(1 + s)]])
display(N)
Explanation: Interpolation functions and matrices
End of explanation
H = zeros(2, 8)
for cont in range(4):
    H[0, 2*cont] = N[cont]
    H[1, 2*cont + 1] = N[cont]
display(H.T)
Explanation: The interpolation matrix is a matrix with the interpolation
functions reordered for each degree of freedom.
End of explanation
dHdr = zeros(2, 4)
for cont in range(4):
    dHdr[0, cont] = diff(N[cont], r)
    dHdr[1, cont] = diff(N[cont], s)
display(dHdr)
Explanation: The local derivatives matrix is formed with the derivatives of the interpolation functions
End of explanation
def gauss_int2d(f, x, y):
    acu = 0
    pts = [-1/sqrt(3), 1/sqrt(3)]
    w = [1, 1]
    for i in range(2):
        for j in range(2):
            acu += f.subs({x: pts[i], y: pts[j]})*w[i]*w[j]
    
    return acu
Explanation: Gauss integration
End of explanation
def jaco(dHdr, coord_el):
    return simplify(dHdr * coord_el)
def jaco_inv(dHdr, coord_el):
    jac = jaco(dHdr, coord_el)
    return Matrix([[jac[1, 1], -jac[0, 1]], [-jac[1, 0], jac[0, 0]]])/jac.det()
Explanation: Local matrices generation
The Jacobian of the transformation is given by the product of the derivative
matrix and the coordinates of each element
End of explanation
def B_mat(dHdr, coord_el):
    dHdx = jaco_inv(dHdr, coord_el) * dHdr
    B = zeros(3, 8)
    for cont in range(4):
        B[0, 2*cont] = dHdx[0, cont]
        B[1, 2*cont + 1] = dHdx[1, cont]
        B[2, 2*cont] = dHdx[1, cont]
        B[2, 2*cont + 1] = dHdx[0, cont]
    return simplify(B)
Explanation: We can re-arrange the derivatives of the interpolation function as a matrix that
translates from displacements to strains, commonly denoted by $B$
End of explanation
def local_mass(H, coord_el, rho):
    det = jaco(dHdr, coord_el).det()
    integrand = rho * det * expand(H.T * H) 
    mass_mat = zeros(8, 8)
    for row in range(8):
        for col in range(row, 8):
            mass_mat[row, col] = gauss_int2d(integrand[row, col], r, s)
            mass_mat[col, row] = mass_mat[row, col]
    return mass_mat
def local_stiff(dHdr, coord_el, C):
    det = jaco(dHdr, coord_el).det()
    B = B_mat(dHdr, coord_el)
    integrand = det * expand(B.T * C * B)
    stiff_mat = zeros(8, 8)
    for row in range(8):
        for col in range(row, 8):
            stiff_mat[row, col] = gauss_int2d(integrand[row, col], r, s)
            stiff_mat[col, row] = stiff_mat[row, col]
    return stiff_mat
Explanation: With these elements we can form the local stiffness and mass matrices.
End of explanation
def assembler(coords, elems, mat_props, const="plane_stress"):
    ncoords = coords.shape[0]
    stiff_glob = zeros(2*ncoords, 2*ncoords)
    mass_glob = zeros(2*ncoords, 2*ncoords)
    for el_cont, elem in enumerate(elems):
        E, nu, rho = mat_props[el_cont]
        C = const_model(E, nu, const=const)
        coord_el = coords[elem, :]
        stiff_loc = local_stiff(dHdr, coord_el, C)
        mass_loc = local_mass(H, coord_el, rho)
        for row in range(4):
            for col in range(4):
                row_glob = elem[row]
                col_glob = elem[col]
                # Stiffness matrix
                stiff_glob[2*row_glob, 2*col_glob] += stiff_loc[2*row, 2*col]
                stiff_glob[2*row_glob, 2*col_glob + 1] += stiff_loc[2*row, 2*col + 1]
                stiff_glob[2*row_glob + 1, 2*col_glob] += stiff_loc[2*row + 1, 2*col]
                stiff_glob[2*row_glob + 1, 2*col_glob + 1] += stiff_loc[2*row + 1, 2*col + 1]
                # Mass matrix
                mass_glob[2*row_glob, 2*col_glob] += mass_loc[2*row, 2*col]
                mass_glob[2*row_glob, 2*col_glob + 1] += mass_loc[2*row, 2*col + 1]
                mass_glob[2*row_glob + 1, 2*col_glob] += mass_loc[2*row + 1, 2*col]
                mass_glob[2*row_glob + 1, 2*col_glob + 1] += mass_loc[2*row + 1, 2*col + 1]
    return stiff_glob, mass_glob
Explanation: Assembly process
End of explanation
coords = Matrix([
    [-1, -1],
    [1, -1],
    [1, 1],
    [-1, 1]])
elems = [[0, 1, 2, 3]]
mat_props = [[S(8)/3, S(1)/3, 1]] 
stiff, mass = assembler(coords, elems, mat_props, const="plane_strain")
stiff
mass
Explanation: Example: Single element
End of explanation
coords = Matrix([
        [-1, -1],
        [0, -1],
        [1, -1],
        [-1,  0],
        [0,  0],
        [1,  0],
        [-1,  1],
        [0,  1],
        [1,  1]])
elems = [[0, 1, 4, 3],
        [1, 2, 5, 4],
        [3, 4, 7, 6],
        [4, 5, 8, 7]]
mat_props = [[16, S(1)/3, 1]]*4
stiff, _ = assembler(coords, elems, mat_props)
stiff_exact = Matrix([
    [8, 3, -5, 0, 0, 0, 1, 0, -4, -3, 0, 0, 0, 0, 0, 0, 0, 0], 
    [3, 8, 0, 1, 0, 0, 0, -5, -3, -4, 0, 0, 0, 0, 0, 0, 0, 0], 
    [-5, 0, 16, 0, -5, 0, -4, 3, 2, 0, -4, -3, 0, 0, 0, 0, 0, 0], 
    [0, 1, 0, 16, 0, 1, 3, -4, 0, -10, -3, -4, 0, 0, 0, 0, 0, 0], 
    [0, 0, -5, 0, 8, -3, 0, 0, -4, 3, 1, 0, 0, 0, 0, 0, 0, 0], 
    [0, 0, 0, 1, -3, 8, 0, 0, 3, -4, 0, -5, 0, 0, 0, 0, 0, 0], 
    [1, 0, -4, 3, 0, 0, 16, 0, -10, 0, 0, 0, 1, 0, -4, -3, 0, 0], 
    [0, -5, 3, -4, 0, 0, 0, 16, 0, 2, 0, 0, 0, -5, -3, -4, 0, 0], 
    [-4, -3, 2, 0, -4, 3, -10, 0, 32, 0, -10, 0, -4, 3, 2, 0, -4, -3], 
    [-3, -4, 0, -10, 3, -4, 0, 2, 0, 32, 0, 2, 3, -4, 0, -10, -3, -4], 
    [0, 0, -4, -3, 1, 0, 0, 0, -10, 0, 16, 0, 0, 0, -4, 3, 1, 0], 
    [0, 0, -3, -4, 0, -5, 0, 0, 0, 2, 0, 16, 0, 0, 3, -4, 0, -5], 
    [0, 0, 0, 0, 0, 0, 1, 0, -4, 3, 0, 0, 8, -3, -5, 0, 0, 0], 
    [0, 0, 0, 0, 0, 0, 0, -5, 3, -4, 0, 0, -3, 8, 0, 1, 0, 0], 
    [0, 0, 0, 0, 0, 0, -4, -3, 2, 0, -4, 3, -5, 0, 16, 0, -5, 0], 
    [0, 0, 0, 0, 0, 0, -3, -4, 0, -10, 3, -4, 0, 1, 0, 16, 0, 1], 
    [0, 0, 0, 0, 0, 0, 0, 0, -4, -3, 1, 0, 0, 0, -5, 0, 8, 3], 
    [0, 0, 0, 0, 0, 0, 0, 0, -3, -4, 0, -5, 0, 0, 0, 1, 3, 8]])
stiff_exact - stiff
Explanation: Example: $2 \times 2$ mesh
End of explanation
coords = Matrix([
        [0, 0],
        [1, 0],
        [2, 0],
        [0, 1],
        [1, 1],
        [2, 1]])
elems = [[0, 1, 4, 3],
        [1, 2, 5, 4]]
mat_props = [[1, S(3)/10, 1]]*4
stiff, mass = assembler(coords, elems, mat_props, const="plane_stress")
load = zeros(12, 1)
load[5] = -S(1)/2
load[11] = -S(1)/2
load
stiff2 = stiff.copy()
stiff2[0, :] = eye(12)[0, :]
stiff2[:, 0] = eye(12)[:, 0]
stiff2[1, :] = eye(12)[1, :]
stiff2[:, 1] = eye(12)[:, 1]
stiff2[6, :] = eye(12)[6, :]
stiff2[:, 6] = eye(12)[:, 6]
stiff2[7, :] = eye(12)[7, :]
stiff2[:, 7] = eye(12)[:, 7]
sol = linsolve((stiff2, load))
sol
from IPython.core.display import HTML
def css_styling():
    styles = open('./styles/custom_barba.css', 'r').read()
    return HTML(styles)
css_styling()
Explanation: Example: $2 \times 1$ cantilever beam
End of explanation
<END_TASK>