The paper presents a general boundary element approach for analysis of 3D cracks in anisotropic thermomagnetoelectroelastic solids. Dual boundary integral equations are derived, which kernels are explicitly written. These equations do not contain volume integrals in the absence of distributed body heat and extended body forces, which is advantageous comparing to the existing approaches. The issues on the boundary element solution of these equations are discussed in details. The efficient numerical evaluation of kernels based on the trapezoid rule is proposed. Modified Kutt's quadrature with Chebyshev nodes is derived for integration of singular and hypersingular integrals. Nonlinear polynomial mappings are adopted for smoothing the integrand at the crack front, which is advantageous for accurate evaluation of field intensity factors. Special shape functions are introduced, which account for a square-root singularity of extended stress and heat flux at the crack front. The issues on numerical determination of field intensity factors are discussed. Several numerical examples are presented, which show the efficiency (low computational time and high precision) of the proposed boundary element formulation.