%% 3D: P1 element matrix tetrahedral element

x = sym('x', [4 1]);
y = sym('y', [4 1]);
z = sym('z', [4 1]);

A = [ones(4,1), x , y, z ];

detA = det(A);
Ainv = inv(A);
T = Ainv*detA;      %  T/detA == Ainv

grad_phi1=T(2:4,1)  %  /detA
grad_phi2=T(2:4,2)  %  /detA
grad_phi3=T(2:4,3)  %  /detA
grad_phi4=T(2:4,4)  %  /detA

%% Laplace 3D
%   - laplace u = rhs
clear
syms x y z
% u(x,y,z) = x^2*sin(pi*y)*cos(pi*z)
u(x,y,z) = (x^2*(2*x - 3))*cos(pi*y)*cos(pi*z)

gu(x,y,z) = simplify(gradient(u));
% rhs = -diff(u,2)
% rhs = -divergence(gradient(u))
rhs = -laplacian(u,[x,y,z]);
rhs(x,y,z) = simplify(rhs)