This is often required in the stiffness matrix evaluation. Thanks to [this][1] thread I found out how to perform integrals over 2D surface. I've tried to evaluate the integral of the single shape function and got 0.3333. Now I want to have the expression which is being integrated to be more complicated. Namely I want to perform this integral: $\displaystyle\int_S {\frac{\partial{H_1(e,n)}}{\partial{x}}}^2 + {\frac{\partial{H_1(e,n)}}{\partial{y}}}^2 + {\frac{\partial{H_1(e,n)}}{\partial{z}}}^2\ ds$ The derivatives of the first shape function with respect to {x,y,z} are not the same as the derivatives with respect to {e,n}. I should use the inversed Jacobian matrix to translate the values of the derivatives with respect to {e,n} at the integration point to the derivatives with respect to {x,y,z}. But the Jacobian matrix isn't square in this case. [1]: https://scicomp.stackexchange.com/questions/3439/how-to-integrate-polynomial-expression-over-3d-4-node-element