%% 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
% \frac{d\;\varphi }{d\;x}=\frac{d\;\varphi }{d\;\xi }\underset{J^{-1} }{\underbrace{\frac{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
% 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
% \underline{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:
% \nabla u=\left\lbrack \begin{array}{c}
% \frac{\partial u}{\partial x}\\
% \frac{\partial u}{\partial y}
% \end{array}\right\rbrack \approx \nabla U=\nabla \left(\varphi \cdot u\right)=\nabla \varphi \cdot u
% where the gradient of the basis functions is this matrix
% \nabla \varphi =\left\lbrack \begin{array}{ccc}
% \frac{\partial \varphi_1 }{\partial x} & \cdots  & \frac{\partial \varphi_n }{\partial x}\\
% \frac{\partial \varphi_1 }{\partial y} & \cdots  & \frac{\partial \varphi_n }{\partial y}
% \end{array}\right\rbrack
%
%% Iso-parametric map in 2D
% 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 chainrule
% \left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial x}\\
% \frac{\partial \varphi_i }{\partial y}
% \end{array}\right\rbrack =\left\lbrack \begin{array}{cc}
% \frac{\partial \xi }{\partial x} & \frac{\partial \eta }{\partial x}\\
% \frac{\partial \xi }{\partial y} & \frac{\partial \eta }{\partial y}
% \end{array}\right\rbrack\left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial \xi }\\
% \frac{\partial \varphi_i }{\partial \eta }
% \end{array}\right\rbrack
% It's straighforward 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
% \left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial \xi }\\
% \frac{\partial \varphi_i }{\partial \eta }
% \end{array}\right\rbrack =\underset{\mathit{\mathbf{J}}}{\underbrace{\left\lbrack \begin{array}{cc}
% \frac{\partial x}{\partial \xi } & \frac{\partial x}{\partial \eta }\\
% \frac{\partial y}{\partial \xi } & \frac{\partial y}{\partial \eta }
% \end{array}\right\rbrack } } \left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial x}\\
% \frac{\partial \varphi_i }{\partial y}
% \end{array}\right\rbrack
% So we can get
% \left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial x}\\
% \frac{\partial \varphi_i }{\partial y}
% \end{array}\right\rbrack ={\mathit{\mathbf{J}}}^{-1}\left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial \xi }\\
% \frac{\partial \varphi_i }{\partial \eta }
% \end{array}\right\rbrack or \mathit{\mathbf{B}}={\mathit{\mathbf{J}}}^{-1} \hat{\mathit{\mathbf{B}}}
% \hat{\mathit{\mathbf{B}}} =\left\lbrack \begin{array}{ccc}
% \frac{\partial \varphi_1 }{\partial \xi } & \cdots  & \frac{\partial \varphi_n }{\partial \xi }\\
% \frac{\partial \varphi_1 }{\partial \eta } & \cdots  & \frac{\partial \varphi_n }{\partial \eta }
% \end{array}\right\rbrack
% \mathit{\mathbf{J}}=\left\lbrack \begin{array}{ccc}
% \frac{\partial \varphi_1 }{\partial \xi } & \cdots  & \frac{\partial \varphi_n }{\partial \xi }\\
% \frac{\partial \varphi_1 }{\partial \eta } & \cdots  & \frac{\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
% \nabla 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}
% \frac{\partial U}{\partial x}\\
% \frac{\partial U}{\partial y}
% \end{array}\right\rbrack
% =\nabla {\mathit{\mathbf{U}}\left(\xi ,\eta \right)}_{\;} ={\mathit{\mathbf{B}}}_{\;} \left(\xi ,\eta \right)\;{\mathit{\mathbf{U}}}_{\;} =\left\lbrack \begin{array}{c}
% \frac{\partial U}{\partial x}\left(\xi ,\eta \right)\\
% \frac{\partial U}{\partial y}\left(\xi ,\eta \right)
% \end{array}\right\rbrack
%% The Jacobian Determinant
% From multivariable calculus (https://www.khanacademy.org/math/multivariable-calculus/multivariable-derivatives/jacobian/v/the-jacobian-determinant) 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
% d\;x\;d\;y=\det \;\mathit{\mathbf{J}}\;d\xi d\eta
% so that we avoid integrating over the \mathit{\mathbf{x}}-coordinates
% \int_K f\left(\mathit{\mathbf{x}}\right)\;d\;\mathit{\mathbf{x}}=\int_{\hat{K} } \hat{f} \left(\xi \right)\det \;\mathit{\mathbf{J}}\;d\xi
% were \hat{K} is the parametric element and \hat{f} \left(\xi \right):=f\left(\mathit{\mathbf{x}}\left(\xi \right)\right).
%
%% 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

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

eq =  [0.5,0.8] == xy;
eq.'
[xi1, eta1] = solve(eq,[xi, eta])

% Now we can compute the gradient. We'll start with the Jacobian matrix.
%  \mathit{\mathbf{J}}=\hat{\mathit{\mathbf{B}}} \left(\xi ,\eta \right)\;\mathit{\mathbf{P}}

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

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

J = Bh*P

% Then we can compute the gradient of the basis functions
% \left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial x}\\
% \frac{\partial \varphi_i }{\partial y}
% \end{array}\right\rbrack =\mathit{\mathbf{B}}={\mathit{\mathbf{J}}}^{-1} \hat{\mathit{\mathbf{B}}}

B = J\Bh

% And finally the gradient of the solution field
% \nabla u\overset{\;}{\approx} \nabla {\mathit{\mathbf{U}}}_{\;} ={\mathit{\mathbf{B}}}_{\;} \;{\mathit{\mathbf{U}}}_{\;}

gradU = double(B*u)

% 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 - Quadliteral
% On a quadliteral 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\rbrackgiven 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)
xy = simplify(phi*P).'

eq =  [3;4] == xy

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

% 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)
eta1 = eta1(2)

% Now we can compute the gradient. We'll start with the Jacobian matrix.
%  \mathit{\mathbf{J}}=\hat{\mathit{\mathbf{B}}} \left(\xi ,\eta \right)\;\mathit{\mathbf{P}}


Bh(xi,eta) = [diff(phi,xi)
              diff(phi,eta)]
J = Bh*P

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

Bh = matlabFunction(Bh)
J = Bh(xi1,eta1)*P

% Then we can compute the gradient of the basis functions
% \left\lbrack \begin{array}{c}
% \frac{\partial \varphi_i }{\partial x}\\
% \frac{\partial \varphi_i }{\partial y}
% \end{array}\right\rbrack =\mathit{\mathbf{B}}={\mathit{\mathbf{J}}}^{-1} \hat{\mathit{\mathbf{B}}}

B = J\Bh(xi1,eta1)

% And finally the gradient of the solution field
% \nabla u\overset{\;}{\approx} \nabla {\mathit{\mathbf{U}}}_{\;} ={\mathit{\mathbf{B}}}_{\;} \;{\mathit{\mathbf{U}}}_{\;}

gradU = double(B*u)

% 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)
[GP,GW] = GaussQuadliteral(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

% 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 isocurves, shows the non-linear nature of the interpolant.
% The gradient of the solution is then computed {20}^2 times and ploted to further show the curved effect.
%
%
%

function [GP,GW] = GaussQuadliteral(order)
%GaussQuadliteral Gauss integration scheme for quadliteral 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