Download this page as a Matlab LiveScript

Table of contents

Iso-Parametric Map
Iso-parametric map in 2D
1. The Jacobian Determinant
Example 1 - Triangles

Iso-Parametric Map

Let's say we want to compute the gradient of some function $u$

Derivative $\frac{d\;u}{d\;x}\approx \frac{d\;U}{d\;x}=\frac{d}{d\;x}\left(\varphi \cdot u\right)=\frac{d\;\varphi }{d\;x}\cdot u$

so what is $\frac{d\;\varphi }{d\;x}$ ?

Using the chain rule we can derive

\dfrac{d\;\varphi }{d\;x}=\dfrac{d\;\varphi }{d\;\xi }\underset{J^{-1} }{\underbrace{\dfrac{d\;\xi }{d\;x} } }

From which we identify the Jacobian.

Now the multidimensional case

Let's say we have a 2D element with $n$ coordinates

\mathbf P=\left\lbrack \begin{array}{cc} x_1 & y_1 \\ \vdots & \vdots \\ x_n & y_n \end{array}\right\rbrack

and we have a corresponding solution field

\mathbf u =\left\lbrack \begin{array}{c} u_1 \\ \vdots \\ u_n \end{array}\right\rbrack

We now want to compute the partial, spatial derivatives of this field, i.e., the gradient:

\def\arraystretch{2.5} \nabla \bm u=\left\lbrack \begin{array}{c} \dfrac{\partial u}{\partial x}\\ \dfrac{\partial u}{\partial y} \end{array}\right\rbrack \approx \nabla \bm U= \nabla \left(\bm \varphi \cdot \bm u\right) = \nabla \bm \varphi \cdot \bm u

where the gradient of the basis functions is this matrix

\def\arraystretch{2.5} \nabla \bm \varphi =\left\lbrack \begin{array}{ccc} \dfrac{\partial \varphi_1 }{\partial x} & \cdots & \dfrac{\partial \varphi_n }{\partial x}\\ \dfrac{\partial \varphi_1 }{\partial y} & \cdots & \dfrac{\partial \varphi_n }{\partial y} \end{array}\right\rbrack

Iso-parametric map in 2D

We shall introduce the concept using the linear triangle element as an example, starting with a map from the spatial space to the parametric.

The coordinates of a linear triangle in x- and y-coordinates can be interpolated using the linear basis function by a linear combination:

\bm x \left(\xi,\eta \right)=\sum_i \varphi_i \left(\xi ,\eta \right) \bm x_i

with

\bm{\varphi}=\begin{bmatrix}1-\xi-\eta & \xi & \eta\end{bmatrix}

thus we have

\begin{align*} \bm{x}(\xi,\eta) & =\sum\varphi_{i}(\xi,\eta)\bm{x}_{i}\\ & =(1-\xi-\eta)\bm{x}_{1}+\xi\bm{x}_{2}+\eta\bm{x}_{3}\\ & =\bm{x}_{1}+\xi(\bm{x}_{2}-\bm{x}_{1})+\eta(\bm{x}_{3}-\bm{x}_{1}) \end{align*}

Now in order to compute this gradient we need to utilize the Iso-Parametric map, compute the gradient of the basis function with respect to the parametric coordinates and map these to the spatial coordinates using the Jacobian.

We begin be deriving this relation, using the chain rule

\def\arraystretch{2.5} \left\lbrack \begin{array}{c} \dfrac{\partial \varphi_i }{\partial x}\\ \dfrac{\partial \varphi_i }{\partial y} \end{array}\right\rbrack =\left\lbrack \begin{array}{cc} \dfrac{\partial \xi }{\partial x} & \dfrac{\partial \eta }{\partial x}\\ \dfrac{\partial \xi }{\partial y} & \dfrac{\partial \eta }{\partial y} \end{array}\right\rbrack \left\lbrack \begin{array}{c} \dfrac{\partial \varphi_i }{\partial \xi }\\ \dfrac{\partial \varphi_i }{\partial \eta } \end{array}\right\rbrack

It's straightforward to compute the $\xi$ - and $\eta$ -derivatives but we do not have access to the partial derivatives of $\xi$ with respect to $x$ , but vice versa. So using the chain rule we can compute

\def\arraystretch{2.5} \left\lbrack \begin{array}{c} \dfrac{\partial \varphi_i }{\partial \xi }\\ \dfrac{\partial \varphi_i }{\partial \eta } \end{array}\right\rbrack =\underset{\mathit{\mathbf{J} } }{\underbrace{\left\lbrack \begin{array}{cc} \dfrac{\partial x}{\partial \xi } & \dfrac{\partial y}{\partial \xi }\\ \dfrac{\partial x}{\partial \eta } & \dfrac{\partial y}{\partial \eta } \end{array}\right\rbrack } } \left\lbrack \begin{array}{c} \dfrac{\partial \varphi_i }{\partial x}\\ \dfrac{\partial \varphi_i }{\partial y} \end{array}\right\rbrack

So we can get

\def\arraystretch{2.5} \begin{bmatrix} \dfrac{\partial \varphi_i }{\partial x}\\ \dfrac{\partial \varphi_i }{\partial y} \end{bmatrix} = \bm J^{-1} \begin{bmatrix} \dfrac{\partial \varphi_i }{\partial \xi }\\ \dfrac{\partial \varphi_i }{\partial \eta } \end{bmatrix} \text{or } \mathbf B = \bm J ^{-1} \hat{\mathbf B}

\def\arraystretch{2.5} \hat{\mathbf B } =\left\lbrack \begin{array}{ccc} \dfrac{\partial \varphi_1 }{\partial \xi } & \cdots & \dfrac{\partial \varphi_n }{\partial \xi }\\ \dfrac{\partial \varphi_1 }{\partial \eta } & \cdots & \dfrac{\partial \varphi_n }{\partial \eta } \end{array}\right\rbrack

\def\arraystretch{2.5} \boldsymbol J =\left\lbrack \begin{array}{ccc} \dfrac{\partial \varphi_1 }{\partial \xi } & \cdots & \dfrac{\partial \varphi_n }{\partial \xi }\\ \dfrac{\partial \varphi_1 }{\partial \eta } & \cdots & \dfrac{\partial \varphi_n }{\partial \eta } \end{array}\right\rbrack \left\lbrack \begin{array}{cc} x_1 & y_1 \\ \vdots & \vdots \\ x_n & y_n \end{array}\right\rbrack =\hat{\mathit{\mathbf{B} } } \;\mathit{\mathbf{P} }

Thus our gradient of the solution field is approximated using the derivatives of our basis functions

\def\arraystretch{2.5} \nabla \bm u\overset{2D}{\approx} \nabla {\mathit{\mathbf{U} } }_{2\times 1} ={\mathit{\mathbf{B} } }_{2\times n} \;{\mathit{\mathbf{U} } }_{n\times 1} =\left\lbrack \begin{array}{c} \dfrac{\partial U}{\partial x}\\ \dfrac{\partial U}{\partial y} \end{array}\right\rbrack

\def\arraystretch{2.5} =\nabla {\mathit{\mathbf{U} }\left(\xi ,\eta \right)}_{\;} ={\mathit{\mathbf{B} } }_{\;} \left(\xi ,\eta \right)\;{\mathit{\mathbf{U} } }_{\;} =\left\lbrack \begin{array}{c} \dfrac{\partial U}{\partial x}\left(\xi ,\eta \right)\\ \dfrac{\partial U}{\partial y}\left(\xi ,\eta \right) \end{array}\right\rbrack

Now, for the linear triangle example we get specifically (using the map in (8))

\def\arraystretch{2} \bm J =\begin{bmatrix}\frac{\partial x}{\partial\xi} & \frac{\partial x}{\partial\eta}\\ \frac{\partial y}{\partial\xi} & \frac{\partial y}{\partial\eta} \end{bmatrix}=\begin{bmatrix}x_{2}-x_{1} & x_{3}-x_{1}\\ y_{2}-y_{1} & y_{3}-y_{1} \end{bmatrix}

or using

\hat{\mathbf{B}}=\begin{bmatrix}-1 & 1 & 0\\ -1 & 0 & 1 \end{bmatrix}

\def\arraystretch{1.5} \bm J =\hat{\mathbf{B}}\mathbf{P}=\begin{bmatrix}-1 & 1 & 0\\ -1 & 0 & 1 \end{bmatrix}\begin{bmatrix}x_{1} & y_{1}\\ x_{2} & y_{2}\\ x_{3} & y_{3} \end{bmatrix}=\begin{bmatrix}x_{2}-x_{1} & x_{3}-x_{1}\\ y_{2}-y_{1} & y_{3}-y_{1} \end{bmatrix}

The Jacobian Determinant

From multivariable calculus we know that the determinant of the Jacobian scales the area of the transformation.

We use the determinant of the Jacobian to scale the reference element area to the physical element by

dxdy = \det \bm J d\xi d\eta

so that we avoid integrating over the $x-,y$ -coordinates

\int_K f\left(x,y\right) dx dy = \int_{\hat{K} } \hat{f} \left(\xi, \eta \right) \det \bm J d\xi d\eta

were $\hat{K}$ is the parametric element and $\hat{f} \left(\xi, \eta \right) := f(x(\xi,\eta), y(\xi, \eta))$ .

Example 1 - Triangles

On a triangle element with nodes given by $\mathit{\mathbf{P} }\;$ below, we want to compute the gradient of $u$ in the point $\left(x,y\right)=\left\lbrack 0\ldotp 5,0\ldotp 8\right\rbrack$ given that the nodal values of $u$ are given below. Plot the directions of the gradient vector at that point.

P = [-1 0
     2 0
     1, 2]

u = [-1; 2; 3]

figure
patch(P(:,1),P(:,2),'w','FaceColor','interp','CData',u); hold on
plot(P(:,1),P(:,2),'ok','MarkerFaceColor','k','MarkerSize',8,'DisplayName','Nodes')
xlabel('x')
ylabel('y')
np = size(P,1);
xm = mean(P,1);
Tp = P+(xm-P)*0.2;
text(Tp(:,1),Tp(:,2),cellstr(num2str((1:np)')),'FontSize',14)
colorbar
plot(0.5,0.8,'*b')

Let's compute the corresponding parametric point by solving

\begin{bmatrix}0.5\\ 0.8 \end{bmatrix}=\bm{\varphi}(\xi,\eta)\mathbf{P}

syms xi eta phi
phi(xi,eta) = [1-xi-eta, xi, eta]

\left(\begin{array}{ccc} 1-\xi -\eta & \xi & \eta \end{array}\right)

xy = phi*P

\left(\begin{array}{cc} 2\,\eta +3\,\xi -1 & 2\,\eta \end{array}\right)

eq =  [0.5,0.8] == xy;
eq.'

\left(\begin{array}{c} \frac{1}{2}=2\,\eta +3\,\xi -1\\ \frac{4}{5}=2\,\eta \end{array}\right)

[xi1, eta1] = solve(eq,[xi, eta])

\xi_1 = \frac{7}{30} \text{, } \eta_1=\frac{2}{5}

Now we can compute the gradient. We'll start with the Jacobian matrix.

\bm J =\hat{\mathit{\mathbf{B} } } \left(\xi ,\eta \right)\;\mathit{\mathbf{P} }

Bh = [diff(phi,xi)
      diff(phi,eta)]

\left(\begin{array}{ccc} -1 & 1 & 0\\ -1 & 0 & 1 \end{array}\right)

We see that this is a constant matrix, independent of where we evaluate it in $\left\lbrack \xi ,\eta \right\rbrack$ .

J = Bh*P

\left(\begin{array}{cc} 3 & 0\\ 2 & 2 \end{array}\right)

Then we can compute the gradient of the basis functions

\left\lbrack \begin{array}{c} \dfrac{\partial \varphi_i }{\partial x}\\ \dfrac{\partial \varphi_i }{\partial y} \end{array}\right\rbrack =\mathit{\mathbf{B} }=\bm J^{-1} \hat{\mathit{\mathbf{B} } }

B = J\Bh

\left(\begin{array}{ccc} -\frac{1}{3} & \frac{1}{3} & 0\\ -\frac{1}{6} & -\frac{1}{3} & \frac{1}{2} \end{array}\right)

And finally the gradient of the solution field

\nabla \bm u \approx \nabla {\mathit{\mathbf{U} } }_{\;} ={\mathit{\mathbf{B} } }_{\;} \;{\mathit{\mathbf{U} } }_{\;}

gradU = double(B*u)

gradU = 2x1    
     1
     1

Let's plot the direction

quiver(0.5,0.8,gradU(1),gradU(2),0.8,'k')

Note that no matter which point in the triangle we look at the gradient is the same everywhere. This is always the case with linear triangles. The gradient is always constant.

Example 2 - Quadrilateral

On a quadrilateral element with nodes given by $\mathit{\mathbf{P} }\;$ below, we want to compute the gradient of $u$ in the point $\left(x,y\right)=\left\lbrack 3,4\right\rbrack$ given that the nodal values of $u$ are given below. Plot the directions of the gradient vector at that point.

clear;
P = [-1 -2 
     6 -2
     7 8
     -2 6];
u = [5
     50
     40
     28]

figure
patch(P(:,1),P(:,2),'w','FaceColor','interp','CData',u); hold on
plot(P(:,1),P(:,2),'ok','MarkerFaceColor','k','MarkerSize',8,'DisplayName','Nodes')
xlabel('x')
ylabel('y')
np = size(P,1);
xm = mean(P,1);
Tp = P+(xm-P)*0.2;
text(Tp(:,1),Tp(:,2),cellstr(num2str((1:np)')),'FontSize',14)
axis equal

Let's compute the corresponding parametric point

phi = @(xi,eta)[(eta-1.0).*(xi-1.0),-xi.*(eta-1.0),eta.*xi,-eta.*(xi-1.0)]

syms xi eta
phi = sym(phi)

\left(\begin{array}{cccc} {\left(\eta -1\right)}\,{\left(\xi -1\right)} & -\xi \,{\left(\eta -1\right)} & \eta \,\xi & -\eta \,{\left(\xi -1\right)} \end{array}\right)

xy = simplify(phi*P).'

\left(\begin{array}{c} 7\,\xi -\eta +2\,\eta \,\xi -1\\ 8\,\eta +2\,\eta \,\xi -2 \end{array}\right)

eq =  [3;4] == xy

\left(\begin{array}{c} 3=7\,\xi -\eta +2\,\eta \,\xi -1\\ 4=8\,\eta +2\,\eta \,\xi -2 \end{array}\right)

[xi1, eta1] = solve(eq,[xi, eta]);
xi1 = double(xi1)

xi1 = 2x1    
   -4.8458
    0.5601

eta1 = double(eta1)

eta1 = 2x1    
   -3.5468
    0.6579

We are getting two solutions, one is false. We know that $\left\lbrack \xi ,\eta \right\rbrack \in \left\lbrack 0,1\right\rbrack$ , se we select the second solution.

xi1 = xi1(2)

xi1 = 0.5601

eta1 = eta1(2)

eta1 = 0.6579

Now we can compute the gradient. We'll start with the Jacobian matrix.

\bm J =\hat{\mathbf B} \left(\xi ,\eta \right) \mathbf P

Bh(xi,eta) = [diff(phi,xi)
              diff(phi,eta)]

\left(\begin{array}{cccc} \eta -1 & 1-\eta & \eta & -\eta \\ \xi -1 & -\xi & \xi & 1-\xi \end{array}\right)

J = Bh*P

\left(\begin{array}{cc} 2\,\eta +7 & 2\,\eta \\ 2\,\xi -1 & 2\,\xi +8 \end{array}\right)

We see that $\hat{\mathit{\mathbf{B} } }$ is not constant and as such the Jacobian is not constant. From now on we will evaluate everything using our parametric point.

Bh = matlabFunction(Bh)

Bh = 
    @(xi,eta)reshape([eta-1.0,xi-1.0,-eta+1.0,-xi,eta,xi,-eta,-xi+1.0],[2,4])

J = Bh(xi1,eta1)*P

J = 2x2    
    8.3158    1.3158
    0.1203    9.1203

Then we can compute the gradient of the basis functions

\def\arraystretch{2.5} \left\lbrack \begin{array}{c} \dfrac{\partial \varphi_i }{\partial x}\\ \dfrac{\partial \varphi_i }{\partial y} \end{array}\right\rbrack =\mathit{\mathbf{B} }= \bm J^{-1} \hat{\mathit{\mathbf{B} } }

B = J\Bh(xi1,eta1)

B = 2x4    
   -0.0336    0.0510    0.0695   -0.0869
   -0.0478   -0.0621    0.0605    0.0494

And finally the gradient of the solution field

\nabla \bm u\overset{\;}{\approx} \nabla {\mathit{\mathbf{U} } }_{\;} ={\mathit{\mathbf{B} } } {\mathit{\mathbf{U} } }_{\;}

gradU = double(B*u)

gradU = 2x1    
    2.7281
    0.4592

Let's plot the direction

quiver(3,4,gradU(1),gradU(2),0.8,'k')

Example 3

Let's compute the gradient of the solution in the Gauss points

phi = matlabFunction(phi)

phi = 
    @(eta,xi)[(eta-1.0).*(xi-1.0),-xi.*(eta-1.0),eta.*xi,-eta.*(xi-1.0)]

[GP,GW] = GaussQuadrilateral(2);

So we're using 4 Gauss points here, for no particular reason.

This means we need to compute the Jacobian and everything in these 4 points and plot it. We will do this using a for loop.

for ig = 1:length(GW)
    xi = GP(ig,1);
    eta = GP(ig,2);
    J=Bh(xi,eta)*P;
    B = J\Bh(xi,eta);
    gradU = B*u
    
    xc = phi(xi,eta)*P;

    plot(xc(1),xc(2),'kx')
    quiver(xc(1),xc(2),gradU(1),gradU(2),0.6,'k')
end

gradU = 2x1    
    4.9952
    2.2452

gradU = 2x1    
    1.8390
    2.0288

gradU = 2x1    
    5.1587
   -0.6270

gradU = 2x1    
    2.2956
   -0.4544

Note how the gradient differs in each point. We are visualizing the solution field using patch and the built-in interpolate function, which is linear. In order to see this non-linear field better. We shall introduce a new way of visualizing a color field on these elements. We can use the basis functions to interpolate the solution field

u\left(\xi ,\eta \right)\approx \varphi \left(\xi ,\eta \right)\cdot \mathbf{u}

N = 100;
xi = linspace(0,1,N);
eta = xi;

[Xi, Eta] = meshgrid(xi,eta);
U = phi(Xi(:),Eta(:))*u;
PP = phi(Xi(:),Eta(:))*P;
X = reshape(PP(:,1),N,N);
Y = reshape(PP(:,2),N,N);
U = reshape(U,N,N);

figure; axis equal tight; hold on
surf(X,Y,Y*0,U,'EdgeColor','none')
colormap jet
colorbar
contour(X,Y,U,5,'EdgeColor','k')

N = 20;
xi = linspace(0,1,N);
eta = xi;
[Xi, Eta] = meshgrid(xi,eta);
Xi = Xi(:); Eta = Eta(:);
PP = phi(Xi(:),Eta(:))*P;

GU = zeros(N^2,2);
for i = 1:N^2
    xi = Xi(i);
    eta = Eta(i);
    J=Bh(xi,eta)*P;
    B = J\Bh(xi,eta);
    gradU = B*u;
    GU(i,:) = gradU';
end

quiver(PP(:,1), PP(:,2), GU(:,1), GU(:,2), 0.6, 'k')

Here we plot the solution variable in ${100}^2$ evenly spaced parametric points and map these points to the real space. This shows that the interpolated solution field is curved. Plotting 5 black iso-curves, shows the non-linear nature of the interpolant.

The gradient of the solution is then computed ${20}^2$ times and plotted to further show the curved effect.

function [GP,GW] = GaussQuadrilateral(order)
%GaussQuadrilateral Gauss integration scheme for quadrilateral 2D element
%   [x,y] in [0,1]

[gp,gw] = gauss(order,0,1);

[GPx,GPy] = meshgrid(gp,gp);
[Gw1,Gw2] = meshgrid(gw,gw);
GW = [Gw1(:), Gw2(:)];
GW = prod(GW,2);
GP = [GPx(:), GPy(:)];

end

function [x, w, A] = gauss(n, a, b)
%------------------------------------------------------------------------------
% gauss.m
%------------------------------------------------------------------------------
%
% Purpose:
%
% Generates abscissas and weigths on I = [ a, b ] (for Gaussian quadrature).
%
%
% Syntax:
%
% [x, w, A] = gauss(n, a, b);
%
%
% Input:
%
% n    integer    Number of quadrature points.
% a    real       Left endpoint of interval.
% b    real       Right endpoint of interval.
%
%
% Output:
%
% x    real       Quadrature points.
% w    real       Weigths.
% A    real       Area of domain.
%------------------------------------------------------------------------------
% 3-term recurrence coefficients:
n = 1:(n - 1);
beta = 1 ./ sqrt(4 - 1 ./ (n .* n));
% Jacobi matrix:
J = diag(beta, 1) + diag(beta, -1); 
% Eigenvalue decomposition:
%
% e-values are used for abscissas, whereas e-vectors determine weights.
%
[V, D] = eig(J);
x = diag(D);
% Size of domain:
A = b - a;
% Change of interval:
%
% Initally we have I0 = [ -1, 1 ]; now one gets I = [ a, b ] instead.
%
% The abscissas are Legendre points.
%
if ~(a == -1 && b == 1)
  x = 0.5 * A * x + 0.5 * (b + a);
end
% Weigths:
w = V(1, :) .* V(1, :);
end