I'm writing a code for solving PDEs through the finite element method. In particular, I'm facing with 3D problems, in which I don't know how to calculate shape functions derivatives on the boundaries (3D triangular elements).

In the 2D triangle, the conversion from local to global shape functions derivatives is done through the following transformation:

$$ \frac{\partial\phi}{\partial \xi} = \frac{\partial \phi}{\partial x}\frac{\partial x}{\partial \xi}+\frac{\partial \phi}{\partial y}\frac{\partial y}{\partial \xi}\\ \frac{\partial\phi}{\partial \eta} = \frac{\partial \phi}{\partial x}\frac{\partial x}{\partial \eta}+\frac{\partial \phi}{\partial y}\frac{\partial y}{\partial \eta} $$


$$ \left[ \begin{array}{c} \frac{\partial \phi}{\partial \xi} \\ \frac{\partial \phi}{\partial \eta} \end{array} \right]= \left[ \begin{array}{cc} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta}\\ \end{array} \right] \left[ \begin{array}{c} \frac{\partial \phi}{\partial x}\\ \frac{\partial \phi}{\partial y} \end{array} \right]= J\left[ \begin{array}{c} \frac{\partial \phi}{\partial x}\\ \frac{\partial \phi}{\partial y} \end{array} \right] $$

The inverse becomes:

$$ \left[ \begin{array}{c} \frac{\partial \phi}{\partial x}\\ \frac{\partial \phi}{\partial y} \end{array} \right] = J^{-1} \left[ \begin{array}{c} \frac{\partial \phi}{\partial \xi} \\ \frac{\partial \phi}{\partial \eta} \end{array} \right] $$

However, adding the $z$ coordinate, the Jacobian becomes a $3 \times 2$ matrix, and its inverse doesn't exist.

How can I solve this problem? How can I pass from $\frac{\partial \phi}{\partial \xi}$ and $\frac{\partial \phi}{\partial \eta}$ to $\frac{\partial \phi}{\partial x}$, $\frac{\partial \phi}{\partial y}$, and $\frac{\partial \phi}{\partial z}$?

How is it done traditionally in FEM codes? If any reference or solution can be shown it would be fantastic.

  • 1
    $\begingroup$ Welcome to Computational Science SE! I've added formulas you've written in MathJax. Please, be sure to check if everything still looks good and use MathJax in your future questions and interactions on Computational Science. $\endgroup$
    – Anton Menshov
    Sep 7, 2020 at 16:08
  • 1
    $\begingroup$ For the boundary your normal 3D shape functions should still be valid. The only difference for the boundary should be doing the surface integration correctly with a surface determinant vs the det J $\endgroup$
    – wwfe
    Sep 8, 2020 at 15:06
  • $\begingroup$ Yes, as you tell for the integral I'm using the surface instead of the jacobian. The problem, however, is not related to the shape functions (their values are the same as in the 2d case) but to their derivatives in the global domain (in the local domain they are the same as in the 2d case). In other words I don't know how to pass from the derivatives in the xi, età domain to their values in the x y z domain. Can you help me ti find a solution? $\endgroup$ Sep 8, 2020 at 19:20
  • $\begingroup$ @wwfe, can you expand your comment into an answer? $\endgroup$
    – nicoguaro
    Sep 9, 2020 at 3:21
  • $\begingroup$ Maybe this is too simplistic, but if you recall that every point on a triangle is on a plane, $n_xx+n_yy+n_zz=0$. Use your 2D results for the third transformation: $∂ϕ/∂z=(∂ϕ/∂x)(∂x/∂z)+(∂ϕ/∂y)(∂y/∂z)$ $\endgroup$
    – Charlie S
    Sep 10, 2020 at 13:37

3 Answers 3


For the case of faceted triangle geometry, the derivatives you're looking for ($\frac{d\phi}{dx}$, $\frac{d\phi}{dy}$, $\frac{d\phi}{dz}$) can actually be found without resorting to calculus (chain rule / jacobian), you can deduce them from purely geometrical considerations.

These derivatives are the cartesian (x,y,z) components of the vector function $\nabla \phi$. Since $\phi$ is a linear function, its gradient $\nabla \phi$ is a constant vector, we just need to establish its (m)agnitude and (d)irection $\nabla \phi = m \hat d$ (which do not vary in space). For sake of example, consider the $\phi_0$ function, which has a value of one at $\vec r_0$ and a value of zero at $\vec r_1$ and $\vec r_2$. Recalling that the gradient of a scalar function is orthogonal to its equipotential surfaces, we deduce that the direction $\hat d$ of $\nabla \phi_0$ is orthogonal to the vector that points along the ($\vec r_1$,$\vec r_2$) edge (because this edge is an equipotential surface, the value of $\phi_0$ is uniformly zero along it). You can form a vector orthogonal to this edge by crossing with the surface normal:

n = cross(r1-r0,r2-r0);
n = n / length(n);
d = cross(n, r2-r1);
d = d / length(d);

enter image description here

This establishes the direction $\hat d$. To establish the magnitude $m$, recall the theorem for line integrals: $\int_{\vec a}^{\vec b} \nabla f(\vec r) \cdot d\vec l = f(\vec b) - f(\vec a)$. That is, when we integrate a gradient along a path, it can be reduced to the difference in potential at the endpoints. This is independent of path. A particularly easy path can be found be dropping a perpendicular from $\vec r_0$ to the ($\vec r_1$,$\vec r_2$) edge, calling that point $\vec a$, then integrate along a straight line towards $\vec b = \vec r_0$. Along this path, $\hat d$ and $d\vec l$ are parallel, reducing the dot product to just scalar multiplication. Furthermore, the endpoint values are also known, $\phi_0(\vec a) = 0$ and $\phi_0(\vec b) = 1$. This lets us write: $\int_{\vec a}^{\vec b} m\cdot dl = 1 - 0$, since $m$ doesn't vary in space, pull it out to obtain $m \int_{\vec a}^{\vec b} dl = 1$, so $m$ is just the reciprocal of the length of that path (the altitude of the triangle).

m = 1 / dot(d,r0-r1)

enter image description here

Now in possession of $\nabla \phi_0 = m \hat d$, you can break it back apart into ($\frac{d\phi_0}{dx}$, $\frac{d\phi_0}{dy}$, $\frac{d\phi_0}{dz}$). The other functions $\phi_1$, $\phi_2$ can be found by permuting indices.

Returning to your original framework, the catch is that our basic intuition for the shape functions, that they have value one at $\vec r_0$ and value zero $\vec r_1$, $\vec r_2$ is insufficient. There are infinitely many functions that satisfy this constraint in $R^3$ (basically you can "rotate" or "shear" the whole coordinate system around the $1,2$ axis and still satisfy those constraints). This leads to a mismatch in dimension / singularity in the Jacobian. Another constraint is necessary to fully specify the behavior of $\phi$ in $R^3$, here we have (somewhat tacitly) chosen that $\frac{d\phi}{dn} = 0$, ie the problem remains invariant as you slide up/down the surface normal. As an alternative, you can also write this constraint as an additional row of the Jacobian to make it invertible (once you advance past faceted geometry, towards FEM on curvilinear/NURBS/etc geometry, the simple vector arithmetic approach outlined by this answer is no longer applicable).


The derivatives of the basis functions with respect to the physical space are given by \begin{equation} \phi_{,x}=\phi_{,\xi}\xi_{,x}+\phi_{,\eta}\eta_{,x} \\ \phi_{,y}=\phi_{,\xi}\xi_{,y}+\phi_{,\eta}\eta_{,y} \\ \phi_{,z}=\phi_{,\xi}\xi_{,z}+\phi_{,\eta}\eta_{,z} \end{equation} Hence, you need $\xi_{,x}$, $\xi_{,y}$, $\xi_{,z}$, $\eta_{,x}$, $\eta_{,y}$, $\eta_{,z}$. They can be obtained as follows. First, we write \begin{equation} \begin{bmatrix} dx \\ dy \\ dz \end{bmatrix} = \underbrace{\begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \\ \frac{\partial z}{\partial \xi} & \frac{\partial z}{\partial \eta} \end{bmatrix}}_{\mbox{$L$}} \begin{bmatrix} d\xi \\ d\eta \end{bmatrix} \end{equation} Now we have \begin{equation} L^T \begin{bmatrix} dx \\ dy \\ dz \end{bmatrix} = L^T L \begin{bmatrix} d\xi \\ d\eta \end{bmatrix} \\ ( L^T L)^{-1}L^T \begin{bmatrix} dx \\ dy \\ dz \end{bmatrix} = \begin{bmatrix} d\xi \\ d\eta \end{bmatrix} \end{equation} This can be rewritten with the 2x3 matrix \begin{equation} A = ( L^T L)^{-1}L^T \end{equation} to the following form \begin{equation} A_{11} dx + A_{12} dy + A_{13} dz = d\xi \\ A_{21} dx + A_{22} dy + A_{23} dz = d\eta \end{equation} or, equivalently \begin{equation} \frac{\partial \xi}{\partial x} dx + \frac{\partial \xi}{\partial y} dy + \frac{\partial \xi}{\partial z} dz = d\xi \\ \frac{\partial \eta}{\partial x} dx + \frac{\partial \eta}{\partial y} dy + \frac{\partial \eta}{\partial z} dz = d\eta \end{equation} So finally, we get by comparison \begin{equation} \begin{bmatrix} \frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} & \frac{\partial \xi}{\partial z} \\ \frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} & \frac{\partial \xi}{\partial z} \end{bmatrix} = \begin{bmatrix} \xi_{,x} & \xi_{,y} & \xi_{,z} \\ \eta_{,x} & \eta_{,y} & \eta_{,z} \end{bmatrix} = \begin{bmatrix} A_{11} & A_{12} & A_{13} \\ A_{21}& A_{22} & A_{23} \end{bmatrix} = ( L^T L)^{-1}L^T \end{equation}

and you are able to compute $\phi_{,x}$, $\phi_{,y}$ and $\phi_{,z}$.

I could find this in Section 3.4 of

Gao, X.-W. (1999). 3D non-linear and multi-region boundary element stress analysis. Ph.D. thesis. University of Glasgow.


An easy way to obtain the derivatives is the following approach.

You start with a set of quadrature points in the reference facet (in this case, a reference triangle). Next you map these points to the global domain, onto the correct global facet, using the local-to-global mapping of the corresponding facet.

At this point you must know which element this global facet belongs to and use the inverse local-to-global mapping of the corresponding tetrahedron to map these points to the reference tetrahedron. Now you are in the position to evaluate the shape function derivatives at the boundary of the reference tetrahedron.

Finally, you transform these values back to the global coordinate system. Due to the construction, the values are now related to the correct global facet.

  • $\begingroup$ Thank you for the answer. However, I asked how to calculate the global shape functions derivatives in the case of a 3d triangle, and not of a tethraedron (in this case it is vert simple to shift from the local to global derivatives because the jacobian matrix is a 3vs3 matrix, while in the case of a 3d triangle the jacobian becomes a 3 vs 2 matrix and this simple procedure cannot be applied). Thank you in advance for the answer $\endgroup$ Sep 9, 2020 at 18:13
  • $\begingroup$ Maybe you realize that the approach is probably not correct just by considering the impossibility of inverting $J$. You cannot magically turn two local derivatives into three global derivatives without going through the corresponding tetrahedron (for which the facet is the triangle that you were originally talking about). $\endgroup$
    – knl
    Sep 9, 2020 at 21:07
  • $\begingroup$ Consider, e.g., linear elements. On a triangle you have 3 DOF's. No way can these uniquely define three derivatives. You need the fourth DOF, in this case the fourth interior node of the tetrahedron. $\endgroup$
    – knl
    Sep 9, 2020 at 21:15
  • $\begingroup$ @knl What about triangular elements defined over a surface? Like shell elements of triangles for the Laplace-Beltrami operator over a surface? $\endgroup$
    – nicoguaro
    Sep 10, 2020 at 12:01
  • 1
    $\begingroup$ Sure, but fixing $\frac{\partial \phi}{\partial n}= 0$ such as the answer suggests is not suitable in most scenarios where I've encountered the need of full gradient in boundary forms. Those are: 1) implementing mortar methods, 2) implementing Nitsche's method, 3) dealing with fixed total current BC's in EIT, 4) integrating the gradient over the boundary when evaluating a mesh-dependent norm. I'm not sure what is attempted here but I'm happy to learn. $\endgroup$
    – knl
    Sep 11, 2020 at 19:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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