The $J_k$ integral method is presented to decompose mixed mode stress intensity factors into mode I and mode II component for the anisotropic body with a crack. Narrow rectangular contours introduced by Cherepanov is adapted to obtain a simple form of the $J_k$ integral formula. Substituting functions of a complex variable defining the singular solution near the crack tip in an anisotropic body into the formula, the complete relations between $J_1$, $J_2$ and $K_Ⅰ$, $K_{Ⅱ}$ are obtained. The relation between $J_1$ and $K_Ⅰ$, $K_{Ⅱ}$ is identical to the existing one. To extract mixed mode stress intensity factors from a finite element solution, the $J_k$ integral is evaluated from it. Path independence of the $J_1$ integral is maintained numerically well.
Unlike other path independent integral method developed for mixed mode crack problems the present method does not need calculating the auxiliary solution or decomposing the displacement. Numerical results of problems with existing solutions are compared and demonstrate the accuracy, stability and versatility of the method. The $\hat{J}_k$ integral as an extension of the $J_k$ integral is also presented for the problems with body forces and thermal loads, etc. Numerical results for an arc crack in a rotating disk are given and agree well with the existing Boundary Element results. Numerical results for a cracked plate subjected to thermal loads are also given and agree well with the existing $\hat{J}$ integral results obtained along a contour far from the crack tip and without neglecting the contribution of an area integral.