# Evaluate the local mass matrix over a quadratic curvilinear tetrahedral element

I am evaluating the local mass matrix over a 3D curvilinear 10-node tetrahedral element. According to an index-naming rule, the shape functions are: $$$$\begin{bmatrix} 2L^2_0 - L_0\\ 2L^2_1 - L_1\\ 2L^2_2 - L_2\\ 2L^2_3 - L_3\\ 4L_0L_1\\ 4L_1L_2\\ 4L_0L_2\\ 4L_0L_3\\ 4L_1L_3\\ 4L_2L_3 \end{bmatrix}$$$$. With these shape functions, I am evaluating the local mass matrix over an arbitrary curvilinear tetrahedral element. The mass matrix is defined as: $$$$A_{mn} = \int \psi_m(x,y,z) \psi_n(x,y,z) d\Omega_e = \int \hat{\psi}_m(L_0,L_1,L2)\hat{\psi}_n(L_0,L_1,L2)\ \vert J \vert\ dL_0 dL_1 dL_2$$$$ where $$\Omega_e$$ is the domain of the element, $$\vert J \vert$$ is the determinant of the Jacobian of the geometry transformation. Since there are no derivatives applied, the integrand $$$$\hat{\psi}_m(L_0,L_1,L2)\hat{\psi}_n(L_0,L_1,L2)\ \vert J \vert\$$$$ are high-order polynomial of $$L_0,L_1,L_2$$.

Since this is the curvilinear tetrahedron(edges are not straight), the area/volume coordinates are infeasible. Gauss quadrature rules seem to be the only option. Zienkiewicz's book (The Finite Element Method: Its Basis and Fundamentals) gives the integration points and weights as shown in the attached screenshot:

However, the cubic points/weights work on diffusion terms (due to their lower order of polynomials) but failed on the mass matrix. It is said that these integration points/weights are deduced from (Numerical integration over simplexes and cones, 1956), however, it is difficult for me to digest that paper. Additionally, (Gauss-Legendre quadrature formulas over a tetrahedron, by Rathod, et al.) said higher order integration points/weights are not given in published literature.

Since the mass matrix is commonly used, I believe there is a mature solution to this problem.

My question is: How should I evaluate the local matrix over a quadratic curvilinear tetrahedron?

• Take a look at this reference by Felippa: researchgate.net/publication/…. It has a 14 point Gauss rule for tetrahedra that is useful for integrating the tet10 mass matrix. Apr 2, 2023 at 19:51
• You might also consider subdividing your tetrahedron into four hexahedra, which can be integrated readily using tensor products of Gauss-Legendre quadrature with some suitable jacobian (weight) adjustments to stitch it all together. Though the resulting rules use a lot of points, I feel there is something to be said for a (relatively) simple technique that can be dialed up to arbitrary order of accuracy on demand. Apr 3, 2023 at 1:20

## 1 Answer

John Burkardt, an applied mathematician at University of Pittsburgh, has compiled a huge collection of small programs to tabulate various quadrature rules, including the Felippa rules mentioned in the comments (and other rules suitable for high order tetrahedron integration). You can find these works on personal pages at the various universities he's worked, and it also seems that his work has been duplicated on github in various repositories. See:

https://people.sc.fsu.edu/~jburkardt/m_src/tetrahedron_felippa_rule/tetrahedron_felippa_rule.html

https://people.math.sc.edu/Burkardt/m_src/tetrahedron_felippa_rule/tetrahedron_felippa_rule.html

https://github.com/johannesgerer/jburkardt-m/tree/master/tetrahedron_felippa_rule

Often these rules have been implemented in each of fortran/c/matlab so pick whatever is easiest to work with. You can also just compile them into standalone programs and scrape the abscissa/weights out, hard-coding them into your present application.