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.


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$?


2 Answers 2


I've taken a look on these equations.

I know you have a 1D element, but rewriting the equation would eliminate confusion if you want to write elements in 2D and 3D.

Therefore your stiffness matrix should be

\begin{equation} \mathbf{K}_{el} = \sum_{i}^{n_{gauss}} \mathbf{B}^T \mathbf{E} \mathbf{B} w_{\xi,i} det(J). \end{equation}

for the numerical integration with Gauss. Your stiffness matrix formula used $J$ instead of $det(J)$, which is not wrong. The determinant of the Jacobian is only significant for multidimensional cases.

The formula for your Jacobian should be \begin{equation} J = \sum_{i}^{n} \frac{\partial N_{i}}{\partial \xi}x_{i}. \end{equation}

The Jacobian in the one-dimensional case is defined as $J = \left[ \frac{\partial x}{\partial \xi} \right]$ for \begin{equation} \begin{bmatrix} \frac{\partial N}{\partial \xi} \end{bmatrix} = J \begin{bmatrix} \frac{\partial N}{\partial x} \end{bmatrix} \end{equation}

\begin{equation} \begin{bmatrix} \frac{\partial N}{\partial \xi} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} \end{bmatrix} \begin{bmatrix} \frac{\partial N}{\partial x} \end{bmatrix} \end{equation}

The Jacobian in the isoparametric concept approximates the transformation by the derivatives of the shape function. Here $\frac{\partial x}{\partial \xi}$ has to be approximated. Like the x-coordinate in the element would be calculated by \begin{equation} x = \begin{bmatrix} N_1 & N_2 & N_3 & N_4 \end{bmatrix} \begin{bmatrix} \hat{x}_1 \\ \hat{x}_2 \\ \hat{x}_3 \\ \hat{x}_4 \\ \end{bmatrix} \end{equation} $\frac{\partial x}{\partial \xi}$ is approximated by \begin{equation} \frac{\partial x}{\partial \xi} = \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} \begin{bmatrix} \hat{x}_1 \\ \hat{x}_2 \\ \hat{x}_3 \\ \hat{x}_4 \\ \end{bmatrix}. \end{equation}

Important: The definition of the Jacobian in many books covering finite elements is not consistent to the mathematical definition. So be careful when you are doing your transformations.

When you are doing the transformation of the derivatives of the shape function, you should have \begin{equation} \begin{bmatrix} \frac{\partial N_1}{\partial x_1} & \frac{\partial N_2}{\partial x_2} & \frac{\partial N_3}{\partial x_3} & \frac{\partial N_4}{\partial x_4} \end{bmatrix} = \mathbf{J}^{-1} \begin{bmatrix} \frac{\partial N_1}{\partial \xi_1} & \frac{\partial N_2}{\partial \xi_2} & \frac{\partial N_3}{\partial \xi_3} & \frac{\partial N_4}{\partial \xi_4} \end{bmatrix}, \end{equation}

so that seems to be ok, except that you have the matrix transposed.

Your stiffness matrix then is wrong. Note that you have a dyadic product, where in this case you get a $4 \times 4$ matrix out of $1 \times 4$ vectors. That leads us to

\begin{equation} \mathbf{K}_{el} = \sum_{i}^{n_{gauss}} \begin{bmatrix} \frac{\partial N_1}{\partial x_1} \\ \frac{\partial N_2}{\partial x_2} \\ \frac{\partial N_3}{\partial x_3} \\ \frac{\partial N_4}{\partial x_4} \end{bmatrix} \mathbf{E} \begin{bmatrix} \frac{\partial N_1}{\partial x_1} & \frac{\partial N_2}{\partial x_2} & \frac{\partial N_3}{\partial x_3} & \frac{\partial N_4}{\partial x_4} \end{bmatrix} w_{\xi,i} det(\mathbf{J}) \end{equation} with $w_{\xi,i}$ being the weight for the gauss integration of the integration point $\xi$.

Gauss quadrature for $\mathbf{K}$ is now really simple. Choose the Gauss points in the isoparametric space and the corresponding weights (just look on wikipedia) and use the formulas given.

  • $\begingroup$ Thanks you very much for your answer. Can you explain what is $\omega_{\xi,i}$? Also, I got the definition of Jacobian. But, how would you compute Jacobian matrix? $\endgroup$
    – C34nm
    Commented Aug 25, 2016 at 0:47
  • $\begingroup$ $w_{xi,i}$ is the weight of the Gauss quadrature. I've already given the formula to calculate the Jacobian in your case: \begin{equation} J = \sum_{i}^{n} \frac{\partial N_{i}}{\partial \xi_i} x_{i}. \end{equation} I will edit my answer to make it a little bit clearer. $\endgroup$
    – P. G.
    Commented Aug 25, 2016 at 17:46
  • $\begingroup$ Thanks. Just to clarify a notation you used, $ \hat{x_i}$ are physical (global) nodal coordinates. Correct? $\endgroup$
    – C34nm
    Commented Aug 26, 2016 at 0:35
  • $\begingroup$ Right. That way we have a distinct difference between nodal and general coordinates. $\endgroup$
    – P. G.
    Commented Aug 26, 2016 at 7:31
  • $\begingroup$ Please share if you have a good reference for computational aspects of FEM. I searched several texts, but didn't find much that clearly discusses programming aspects. For example, a number of texts discuss 1D model but it becomes confusing to generalize. Many other texts discuss jacobian integral, but it is not clear how to use it in program. $\endgroup$
    – C34nm
    Commented Aug 29, 2016 at 0:15

From $x(\xi) =N(\xi)x_{i}$ the jacobian is given by:

\begin{equation} J = \frac{\partial{x}}{\partial{\xi}} = \Sigma\frac{\partial{N_{i}}}{\partial{\xi}}x_{i} \end{equation}

The integration itself is done at quadrature points (and thus the jacobian is found by evaluating $\frac{\partial{N}}{\partial{\xi}}$ at quadrature points)

A good reference is Introduction to the Finite Element Method.

Here is some more detail using my own code (C++). Consider a shape function for a quadratic 1D line given by:

void LineQuad1D::shape(double xi)
  //shape functions for quadratic 1D lines at node (xi)
  N[0] = -0.5*xi*(1-xi);
  N[1] = 0.5*xi*(1+xi);
  N[2] = 1-xi*xi;

and corresponding derivative:

void LineQuad1D::dshape(double xi)
  //derivatives of shape functions for quadratic 1D lines at node (xi)
  dN[0] = -0.5+xi; dN[1] = 0.5+xi; dN[2] = -2.0*xi;

Now finding the jacobian at point might look like:

double Element::jacobian(double xi)
  Shape *sh = ShapeFactory::NewShape(Type);

  jmat[0][0] = 0.0;
  for(int k=0;k<Npe;k++){
    jmat[0][0] += (*sh).dN[k]*xpts[k];

  //jacobian, jacobian matrix, and inverse jacobian matrix
  jac = jmat[0][0];  //in your case the dimenion is 1, so you wouldn't need a 2D array for the jacobian
  imat[0][0] = 1.0/jac;

Finding the $B$ matrix at a point looks like:

double Element::tempDiffMatrix(double xi)
  Shape *sh = ShapeFactory::NewShape(Type);


  //set jacobian and inverse jacobian matrix and return jacobian
  double jac = jacobian(xi);

  //temperature differentiation matrix B
  for(int i=0;i<Npe;i++){
    bmat[i] = imat[0][0]*(*sh).dN[i];

  return jac;

And finally stiffness matrix:

void Element::stiffnessMatrix()
  for(int i=0;i<Npe;i++){
    for(int j=0;j<Npe;j++){
      kmat[i][j] = 0.0; //stiffness matrix

  double u=1.0

  QuadratureRule g(nGauss,nDim,Type);
  for(int ip=0;ip<g.nIntPoints;ip++) //loop through quadrature points

    double det = tempDiffMatrix(u);  //creates your B matrix and returns Jacobian, both at integration points xii[ip]
    double dv = det*g.wi[ip];

    for(int i=0;i<Npe;i++){   // Npe = "Number of Points per Element
      for(int j=0;j<Npe;j++){
        double s = 0.0;
        for(int k=0;k<nDim;k++)  // nDim = 1 in your case
          s += bmat[k][i]*bmat[k][j];
        kmat[i][j] += dv*s;

My code actually solves for 1, 2, and 3 dimensions but I have tried to edit out the multi-dimensionality for you since it looks like you are only worried about 1D. Hopefully you can parse the basic idea from these code snippets. Also note that the these code snippets are for solving the heat equation (i.e. no $E$ matrix), but again hopefully you can parse the basic idea.


Ok so for more detail on the Gauss quadrature part here is an example found in "Programming Finite Elements in Java" by Gennadiy Nikishkov. Consider trying to integrate the following:

\begin{equation} I = \int_{-1}^{1}f(\xi)d\xi, \end{equation}

using two integration points. Because we are using two integration points this should integrate to a third degree polynomial. We write the integral as:

\begin{equation} I = f(\xi_{1})w_{1} + f(\xi_{2})w_{2} \end{equation}

where $w_{i}$ are the weights and $\xi_{1}$ are the undetermined points. This formula should be valid for any third degree polynomial, in particular it should work for $f(\xi)=1$, $f(\xi)=\xi$, $f(\xi)=\xi^{2}$, and $f(\xi)=\xi^{3}$:

\begin{equation} \int_{-1}^{1}d\xi = 2 = w_{1} + w_{2}, \end{equation}

\begin{equation} \int_{-1}^{1}\xi d\xi = 0 = \xi_{1}w_{1} + \xi_{2}w_{2}, \end{equation}

\begin{equation} \int_{-1}^{1}\xi^{2} d\xi = 2/3 = \xi_{1}^{2}w_{1} + \xi_{2}^{2}w_{2}, \end{equation}

\begin{equation} \int_{-1}^{1}\xi^{3} d\xi = 0 = \xi_{1}^{3}w_{1} + \xi_{2}^{3}w_{2}, \end{equation}

Solving these results in $\xi_{1}=\frac{1}{\sqrt{3}}$, $\xi_{2}=-\frac{1}{\sqrt{3}}$, $w_{1} = w_{2} = 1$. In general for the one dimensional case:

\begin{equation} I = \int_{-1}^{1}f(\xi)d\xi = \Sigma_{i=1}^{n}f(\xi_{i})w_{i}, \end{equation}

where n is the number of integration points, $\xi_{i}$ the undetermined points, and $w_{i}$ the weights. So in my stiffnessMatrix() function/method above you can see that in pseudo code I am doing the following:

for each integration point p:
    find the Jacobian, dv, at undetermined point xi[p]
    find weight corresponding to undetermined point xi[p]
    find det = jacobian*weight corresponding to point xi[p] 
    find B matrix at the undetermined point xi[p]
    for i = 1:number of ponts in element:
        for j = number of points in element:
          K[i][j] += dv*B[i]B[j]

These undetermined points and weights would be stored in an array or table.

  • $\begingroup$ Thanks James for your answer. In your answer, $x(\xi)=N(\xi)x_i$ you mean global coordinates by $xi$ (not in isoparametric system). Also, I edited the question. Can you comment on that? $\endgroup$
    – C34nm
    Commented Aug 19, 2016 at 0:27
  • $\begingroup$ @C34nm Yes $x_{i}$ are the global coordinates. Looking at your $K$ matrix this looks correct. The link I provided is from a book by Gennadiy Nikishkov titled "Programming Finite Elements in Java". Unfortunately that pdf link does not include the chapter on Gauss quadrature (which is really clear and easy to understand). Fortunately I have the book and I will add an example on this subject tomorrow...its late where I am and I have to work in the morning. $\endgroup$
    – James
    Commented Aug 19, 2016 at 6:33
  • $\begingroup$ @C34nm Ok I added more information which hopefully helps. $\endgroup$
    – James
    Commented Aug 20, 2016 at 3:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.