---

Geometric Definitions

Figure 1 illustrates the quadratic triangular element (TRI6), which contains six nodes: three corner (vertex) nodes and three midside nodes located at the midpoints of each edge. The local coordinate system follows the standard counterclockwise convention \((i,j,k)\), ensuring positive elemental area. The geometric coefficients \((b_i,b_j,b_k)\) and \((c_i,c_j,c_k)\) are defined from the nodal coordinates \((x_i,y_i)\) and are used to compute the derivatives of the shape functions. This element provides quadratic interpolation for displacement, velocity, or pressure fields, and is commonly employed for isoparametric formulations in incompressible flow and structural mechanics applications.

The shape functions \(N_i\) of the quadratic triangular element (TRI6) are constructed from the barycentric coordinates \((L_1,L_2,L_3)\), which satisfy \(L_1 + L_2 + L_3 = 1\). Each vertex function (\(N_1\), \(N_2\), \(N_3\)) provides quadratic interpolation along the corresponding nodal direction, while the midside functions (\(N_4\), \(N_5\), \(N_6\)) ensure \(C^0\) continuity along the edges. These shape functions form a complete quadratic basis over the element, allowing accurate representation of both linear and parabolic variations of the field variables within the domain. Their construction guarantees the Kronecker delta property at the element nodes and the partition of unity over the entire element.

\[\begin{equation} \begin{aligned} N_i &= L_i(2L_i - 1),\\ N_j &= L_j(2L_j - 1),\\ N_k &= L_l(2L_k - 1),\\ N_l &= 4L_iL_j,\\ N_m &= 4L_jL_k,\\ N_n &= 4L_kL_i. \end{aligned} \end{equation}\]

Let \((b_i,b_j,b_k)\) and \((c_i,c_j,c_k)\) be the classical geometric coefficients: \[b_i = y_j - y_k, \qquad b_j = y_k - y_i, \qquad b_k = y_i - y_j,\] \[c_i = x_k - x_j, \qquad c_j = x_i - x_k, \qquad c_k = x_j - x_i,\] and the elemental area \[A = \frac{1}{2} \begin{vmatrix} 1 & x_i & y_i \\ 1 & x_j & y_j \\ 1 & x_k & y_k \end{vmatrix}.\]

Mass Matrix

The element mass matrix \(\mathbf{m}^e\) represents the consistent discretization of the inertial term within the finite element formulation. It is obtained by integrating the product of shape functions over the element domain. For the quadratic triangular element (TRI6), this matrix ensures proper energy conservation and accurate transient response by distributing the mass according to the spatial interpolation of the shape functions. The resulting matrix is symmetric, positive definite, and satisfies the partition of unity, guaranteeing consistent mass representation across the finite element mesh.

\[\begin{equation} \mathbf{m}^e = \int_{\Omega^e} N_i N_j d\Omega^e \end{equation}\]

\[\begin{equation} \mathbf{m}^e = \frac{A}{180} \begin{bmatrix} 6 & -1 & -1 & 0 & -4 & 0 \\ -1 & 6 & -1 & 0 & 0 & -4 \\ -1 & -1 & 6 & -4 & 0 & 0 \\ 0 & 0 & -4 & 32 & 16 & 16 \\ -4 & 0 & 0 & 16 & 32 & 16 \\ 0 & -4 & 0 & 16 & 16 & 32 \end{bmatrix}_{(6,6)} \end{equation}\]

Stiffness Matrices

Matrices \(\mathbf{k}_x^e\) and \(\mathbf{k}_y^e\)

The element stiffness matrices \(\mathbf{k}x^e\) and \(\mathbf{k}y^e\) arise from the spatial derivatives of the shape functions and represent the diffusive or elastic contributions in the finite element formulation. They are defined as

\[\begin{equation} \mathbf{k}_x^e = \int_{\Omega^e} \dfrac{\partial N_i}{\partial x} \dfrac{\partial N_j}{\partial x} d\Omega^e \qquad \mathbf{k}_y^e = \int_{\Omega^e} \dfrac{\partial N_i}{\partial y} \dfrac{\partial N_j}{\partial y} d\Omega^e \end{equation}\]

and collectively form the isotropic stiffness operator \(\mathbf{k}^e = \mathbf{k}_x^e + \mathbf{k}_y^e\). For quadratic triangular elements (TRI6), these matrices provide a second-order accurate spatial approximation of gradients within the element. The resulting matrices are symmetric, positive definite, and ensure correct energy consistency between the nodal degrees of freedom.

\[\begin{equation} \scriptsize \mathbf{k}_x^e = \frac{1}{12A} \begin{bmatrix} 3b_i^2 & -b_i b_j & -b_i b_k & 4b_i b_j & 0 & 4b_i b_k \\ -b_i b_j & 3b_j^2 & -b_j b_k & 4b_i b_j & 4b_j b_k & 0 \\ -b_i b_k & -b_j b_k & 3b_k^2 & 0 & 4b_j b_k & 4b_i b_k \\ 4b_i b_j & 4b_i b_j & 0 & 8(b_i^2+b_i b_j+b_j^2) & 4(b_j^2+b_j b_k+b_i b_j+2b_i b_k) & 4(b_i^2+b_i b_k+b_i b_j+2b_j b_k) \\ 0 & 4b_j b_k & 4b_j b_k & 4(b_i b_j+2b_i b_k+b_j^2+b_j b_k) & 8(b_j^2+b_j b_k+b_k^2) & 4(2b_i b_j+b_i b_k+b_j b_k+b_k^2) \\ 4b_i b_k & 0 & 4b_i b_k & 4(b_i^2+b_i b_j+b_i b_k+2b_j b_k) & 4(2b_i b_j+b_i b_k+b_j b_k+b_k^2) & 8(b_i^2+b_i b_k+b_k^2) \end{bmatrix}_{(6,6)} \end{equation}\]

\[\begin{equation} \scriptsize \mathbf{k}_y^e = \frac{1}{12A} \begin{bmatrix} 3c_i^2 & -c_i c_j & -c_i c_k & 4c_i c_j & 0 & 4c_i c_k \\ -c_i c_j & 3c_j^2 & -c_j c_k & 4c_i c_j & 4c_j c_k & 0 \\ -c_i c_k & -c_j c_k & 3c_k^2 & 0 & 4c_j c_k & 4c_i c_k \\ 4c_i c_j & 4c_i c_j & 0 & 8(c_i^2+c_i c_j+c_j^2) & 4(c_j^2+c_j c_k+c_i c_j+2c_i c_k) & 4(c_i^2+c_i c_k+c_i c_j+2c_j c_k) \\ 0 & 4c_j c_k & 4c_j c_k & 4(c_i c_j+2c_i c_k+c_j^2+c_j c_k) & 8(c_j^2+c_j c_k+c_k^2) & 4(2c_i c_j+c_i c_k+c_j c_k+c_k^2) \\ 4c_i c_k & 0 & 4c_i c_k & 4(c_i^2+c_i c_j+c_i c_k+2c_j c_k) & 4(2c_i c_j+c_i c_k+c_j c_k+c_k^2) & 8(c_i^2+c_i c_k+c_k^2) \end{bmatrix}_{(6,6)} \end{equation}\]

Matrix \(\mathbf{k}_{xy}^e\)

The mixed stiffness matrix \(\mathbf{k}_{xy}^e\) accounts for the coupling between derivatives of the shape functions in the \(x\) and \(y\) directions. It is defined as

\[\begin{equation} \mathbf{k}_{xy}^e = \int_{\Omega^e} \dfrac{\partial N_i}{\partial x} \dfrac{\partial N_j}{\partial y} d\Omega^e \end{equation}\]

and represents the cross-diffusion or shear interaction terms that appear in multidimensional problems. This matrix, together with \(\mathbf{k}_x^e\) and \(\mathbf{k}_y^e\), forms the complete stiffness representation of the element. In isotropic cases, \(\mathbf{k}_{xy}^e\) contributes to maintaining the correct balance of directional derivatives and ensures that the numerical discretization properly captures coupling effects between the two spatial components. The resulting matrix is generally nonsymmetric but, when combined with its transpose, contributes symmetrically to the global stiffness operator, preserving overall consistency and energy conservation.

\[\begin{equation} \scriptsize \mathbf{k}_{xy}^e = \frac{1}{12A} \begin{bmatrix} 3b_i c_i & -b_i c_j & -b_i c_k & 4b_i c_j & 0 & 4b_i c_k \\ -b_j c_i & 3b_j c_j & -b_j c_k & 4b_j c_i & 4b_j c_k & 0 \\ -b_k c_i & -b_k c_j & 3b_k c_k & 0 & 4b_k c_j & 4b_k c_i \\ 4b_j c_i & 4b_i c_j & 0 & 4(2b_i c_i + b_i c_j + b_j c_i + 2b_j c_j) & 4(b_i c_j + 2b_i c_k + b_j c_j + b_j c_k) & 4(b_i c_i + b_i c_k + b_j c_i + 2b_j c_k) \\ 0 & 4b_k c_j & 4b_j c_k & 4(b_j c_i + b_j c_j + 2b_k c_i + b_k c_j) & 4(2b_j c_j + b_j c_k + b_k c_j + 2b_k c_k) & 4(2b_j c_i + b_j c_k + b_k c_i + b_k c_k) \\ 4b_k c_i & 0 & 4b_i c_k & 4(b_i c_i + b_i c_j + b_k c_i + 2b_k c_j) & 4(2b_i c_j + b_i c_k + b_k c_j + b_k c_k) & 4(2b_i c_i + b_i c_k + b_k c_i + 2b_k c_k) \end{bmatrix}_{(6,6)} \end{equation}\]

Gradient and Convective Operators

The gradient operators \(\mathbf{g}_x^e\) and \(\mathbf{g}_y^e\) establish the coupling between the scalar and vector fields within the element, and are defined by

\[\begin{equation} \mathbf{g}_{x}^e = \int_{\Omega^e} \dfrac{\partial N_i}{\partial x} L_j d\Omega^e \qquad \mathbf{g}_{y}^e = \int_{\Omega^e} \dfrac{\partial N_i}{\partial y} L_j d\Omega^e \end{equation}\]

where \(N_i\) and \(L_j\) denote, respectively, the linear shape functions of the velocity and pressure spaces. These matrices are fundamental in mixed finite element formulations, as they link the continuity and momentum equations through the pressure–velocity coupling.

\[\begin{equation} \mathbf{g}_x^e = \frac{1}{6} \begin{bmatrix} b_i & 0 & 0 \\ 0 & b_j & 0 \\ 0 & 0 & b_k \\ b_i+2b_j & 2b_i+b_j & b_i+b_j \\ b_j+b_k & b_j+2b_k & 2b_j+b_k \\ b_i+2b_k & b_i+b_k & 2b_i+b_k \end{bmatrix}_{(6,3)}, \qquad \mathbf{g}_y^e = \frac{1}{6} \begin{bmatrix} c_i & 0 & 0 \\ 0 & c_j & 0 \\ 0 & 0 & c_k \\ c_i+2c_j & 2c_i+c_j & c_i+c_j \\ c_j+c_k & c_j+2c_k & 2c_j+c_k \\ c_i+2c_k & c_i+c_k & 2c_i+c_k \end{bmatrix}_{(6,3)} \end{equation}\]

The convective operators \(\mathbf{gv}_x^e\) and \(\mathbf{gv}_y^e\), on the other hand, represent the advective transport of field quantities and are given by

\[\begin{equation} \mathbf{gv}_{x}^e = \int_{\Omega^e} \dfrac{\partial N_i}{\partial x} N_j d\Omega^e \qquad \mathbf{gv}_{y}^e = \int_{\Omega^e} \dfrac{\partial N_i}{\partial y} N_j d\Omega^e \end{equation}\]

These matrices describe the nonlinear convection terms that arise in advection–diffusion and Navier–Stokes equations, ensuring a consistent discretization of the transport phenomena inside each element. Together, the gradient and convective operators provide the essential link between differential and integral formulations, preserving both conservation and stability properties in the finite element solution.

\[\begin{equation} \mathbf{gv}_x^e = \frac{1}{30} \begin{bmatrix} 2b_i & -b_j & -b_k & -b_i+2b_j & -(b_j+b_k) & -b_i+2b_k \\ -b_i & 2b_j & -b_k & 2b_i-b_j & -b_j+2b_k & -(b_i+b_k) \\ -b_i & -b_j & 2b_k & -(b_i+b_j) & 2b_j-b_k & 2b_i-b_k \\ 3b_i & 3b_j & -b_k & 8(b_i+b_j) & 4(b_j+2b_k) & 4(b_i+2b_k) \\ -b_i & 3b_j & 3b_k & 4(2b_i+b_j) & 8(b_j+b_k) & 4(2b_i+b_k) \\ 3b_i & -b_j & 3b_k & 4(b_i+2b_j) & 4(2b_j+b_k) & 8(b_i+b_k) \end{bmatrix}_{(6,6)} \end{equation}\]

\[\begin{equation} \mathbf{gv}_y^e = \frac{1}{30} \begin{bmatrix} 2c_i & -c_j & -c_k & -c_i+2c_j & -(c_j+c_k) & -c_i+2c_k \\ -c_i & 2c_j & -c_k & 2c_i-c_j & -c_j+2c_k & -(c_i+c_k) \\ -c_i & -c_j & 2c_k & -(c_i+c_j) & 2c_j-c_k & 2c_i-c_k \\ 3c_i & 3c_j & -c_k & 8(c_i+c_j) & 4(c_j+2c_k) & 4(c_i+2c_k) \\ -c_i & 3c_j & 3c_k & 4(2c_i+c_j) & 8(c_j+c_k) & 4(2c_i+c_k) \\ 3c_i & -c_j & 3c_k & 4(c_i+2c_j) & 4(2c_j+c_k) & 8(c_i+c_k) \end{bmatrix}_{(6,6)} \end{equation}\]

Divergent Operator

The discrete divergence operator for the quadratic triangular element (TRI6) is obtained from the transpose of the gradient matrices \(\mathbf{g}_x^e\) and \(\mathbf{g}_y^e\). These matrices relate the pressure shape functions to the derivatives of the velocity shape functions, and are essential in the discretization of the incompressible continuity equation.

\[\begin{equation} \mathbf{d}_{x}^e = \int_{\Omega^e} L_i \dfrac{\partial N_j}{\partial x} d\Omega^e \qquad \mathbf{d}_{y}^e = \int_{\Omega^e} L_i \dfrac{\partial N_j}{\partial y} d\Omega^e \end{equation}\]

\[\begin{equation} \mathbf{d}_x^e = \left(\mathbf{g}_x^e\right)^{\mathsf{T}}, \qquad \mathbf{d}_y^e = \left(\mathbf{g}_y^e\right)^{\mathsf{T}} \end{equation}\]

Explicitly,

\[\begin{equation} \small \mathbf{d}_x^e = \frac{1}{6} \begin{bmatrix} b_i & 0 & 0 & b_i+2b_j & b_j+b_k & b_i+2b_k \\ 0 & b_j & 0 & 2b_i+b_j & b_j+2b_k & b_i+b_k \\ 0 & 0 & b_k & b_i+b_j & 2b_j+b_k & 2b_i+b_k \end{bmatrix}_{(3,6)}, \, \mathbf{d}_y^e = \frac{1}{6} \begin{bmatrix} c_i & 0 & 0 & c_i+2c_j & c_j+c_k & c_i+2c_k \\ 0 & c_j & 0 & 2c_i+c_j & c_j+2c_k & c_i+c_k \\ 0 & 0 & c_k & c_i+c_j & 2c_j+c_k & 2c_i+c_k \end{bmatrix}_{(3,6)}. \end{equation}\]

Slip-Gradient Operator

Additionally, the so-called slip gradient operators are introduced to enforce weak boundary conditions or coupling with immersed boundaries. These matrices, \(\mathbf{g_{slip}}_x^e\) and \(\mathbf{g_{slip}}_y^e\), are defined as:

\[\begin{equation} \mathbf{g_{slip}}_{x}^e = \int_{\Omega^e} N_i \dfrac{\partial L_j}{\partial x} d\Omega^e \qquad \mathbf{g_{slip}}_{y}^e = \int_{\Omega^e} N_i \dfrac{\partial L_j}{\partial y} d\Omega^e \end{equation}\]

\[\begin{equation} \mathbf{gx_{slip}^e} = \frac{1}{6} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ b_i & b_j & b_k \\ b_i & b_j & b_k \\ b_i & b_j & b_k \end{bmatrix}_{(6,3)}, \qquad \mathbf{gy_{slip}^e} = \frac{1}{6} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ c_i & c_j & c_k \\ c_i & c_j & c_k \\ c_i & c_j & c_k \end{bmatrix}_{(6,3)}. \end{equation}\]

These operators are particularly useful when modeling velocity-slip or shear-free boundary conditions in finite element formulations for incompressible Navier–Stokes equations.