Finding jacobian matrix of cubic finite element

Question

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

Answer 1

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.

Answer 2

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);
  (*sh).dshape(xi);

  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);

  (*sh).dshape(xi);

  //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
  {
    u=g.xii[ip];

    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.

Edit

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.