Simulation of many practical problems requires to use non-linear formulations with large displacements, large strains and large rotations. It is well known that the use of Trefftz (T-) functions (i. e. the functions satisfying the governing equations inside the domain) as weighting, or interpolation functions leads to more efficient formulations than those obtained by classical methods. In this paper we will show the use of T-functions and especially T-polynomials, Kelvin, or Kupradze and Boussinesq functions (Green functions with singularity points defined outside of the domain) and their combination in connection with the total Lagrangian formulation for multi-domain BEM (reciprocity based FEM) analysis of displacements and for the post-processing phase in the analysis (evaluation of both gradient of displacements and stress fields). The formulation results in non-singular boundary integrals which has numerical advantages over other formulations using singular boundary integral equations.