I am programming cubic element and very confused how to find interpolation function for cubic finite element.
I have element like this in bi-unit isoparametic system :
$-1 \hspace{15pt} -1/3 \hspace{15pt} 1/3 \hspace{20pt} 1$ <- distance
|-------------|------------|-----------|
$1 \hspace{30pt} 2 \hspace{30pt} 3 \hspace{30pt} 4 \hspace{15pt}$ <- Nodes
Physically, the element goes from 1 to 4. Nodes are located at $[1,2,3,4]$. In natural (isoparametric bi-unit system), nodes $[1,2,3,4]$ are located at $[-1,-1/3,1/3,1]$. We have $\textbf{x} = N_1x_1 + N_2x_2 + N_3x_3 + N_4x_4 $. Where, $N_i$ is shape function in geometric coordinates and $x_i$ are nodal locations. Lagrangian interpolation functions for 1D in terms of physical coordinates are: $$N_1(x) = \frac{ (x-x_2) (x-x_3) (x-x_4)} {(x_1-x_2) (x_1-x_3) (x_1-x_4)}$$ $$N_2(x) = \frac{ (x-x_1) (x-x_3) (x-x_4)} {(x_2-x_1) (x_2-x_3) (x_2-x_4)}$$
similarly for other two nodes (I am not writing them).
In terms of natural coordinates,
$$N_1(\xi) = \frac{ (\xi-\xi_2) (\xi-\xi_3) (\xi-\xi_4)} {(\xi_1-\xi_2) (\xi_1-\xi_3) (\xi_1-\xi_4)}$$ $$N_2(\xi) = \frac{ (\xi-\xi_1) (\xi-\xi_3) (\xi-\xi_4)} {(\xi_2-\xi_1) (\xi_2-\xi_3) (\xi_2-\xi_4)}$$
$$\textbf B = \bigg [\frac{\partial N_1}{\partial x}, \frac{\partial N_2}{\partial x}, \frac{\partial N_3}{\partial x}, \frac{\partial N_4}{\partial x} \bigg ] \\ x(\xi) = N(\xi)x_i $$ Jacobian $J$ is given by $$J = \frac{\partial x}{\partial \xi} $$
When we introduce isoparametric coordinates, $$\textbf B = \frac{\partial N_i}{\partial x} = \frac{\partial N_i}{\partial \xi} \frac{\partial \xi}{\partial x} = \frac{\partial N_i}{\partial \xi}\frac{1}{J}$$
Stiffness matrtix is given by $$\textbf K = \int_{-1}^{1} \textbf{ [B]}^{T} \textbf{ E [B] J} d\xi $$
Can someone tell me how to calculate $\textbf B$, $\textbf J$ for above given element? Many references only linear elements, in which we can easily get away with constants. I tried looking some codes, but no luck as many of them still had linear elements. Any references are welcome.
Edit
After James' answer, I modify my question and reflect findings.
Here is $$\frac{\partial N_i}{\partial \xi} = \begin{bmatrix} -27\xi^2/16 + 9\xi/8 + 1/16 & 81\xi^2/16 - 9\xi/8 - 27/16 & -27\xi^2/16 + 9\xi/8 + 1/16 & -27\xi^2/16 + 9\xi/8 + 1/16 & \end {bmatrix} $$
$$x_i= \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end {bmatrix} $$
So, $$J(\xi) = \frac{\partial N_i}{\partial \xi} x_i$$ is function of $\xi$. $$B(\xi) ^{T} = \frac{1}{J} \begin{bmatrix} \frac{\partial N_1}{\partial \xi} & \frac{\partial N_2}{\partial \xi} & \frac{\partial N_3}{\partial \xi} & \frac{\partial N_4}{\partial \xi} & \end{bmatrix} $$
$B$ is also a function of $\xi$.
$$ K = \int_{-1}^{1} \begin{pmatrix} \frac{1}{J} \begin{bmatrix} \frac{\partial N_1}{\partial \xi} & \frac{\partial N_2}{\partial \xi} & \frac{\partial N_3}{\partial \xi} & \frac{\partial N_4}{\partial \xi} & \end{bmatrix} \end{pmatrix} D \begin{pmatrix} \frac{1}{J} \begin{bmatrix} \frac{\partial N_1}{\partial \xi} \\ \frac{\partial N_2}{\partial \xi} \\ \frac{\partial N_3}{\partial \xi} \\ \frac{\partial N_4}{\partial \xi} \\ \end{bmatrix} \end{pmatrix} J d\xi$$
After that, we find $K$ by Gauss integral as Are these formulas correct, especially $K$? Can someone throw light on how to find Gauss quadrature for $K$?