Also, I note that the tangent seems to be missing the ctan(2) to pick up the effect of the time integrator you use.
For the transient heat conduction, ctan(2) is applied in thtrans3d at line 134 (FEAP 8.4).
Also, one cannot compute the r(4,ii) as shown as the s(i1,j1) will have other terms in it from the thermal gradient before the call.
As far as I see, the DOF 4 components in s are still zero, when it comes to thtrans3d. Before calling thtrans3d in fld3d1tm, only mechanical inertia and geometric tangent terms are added to DOF 1-3 of s.
I have created a 3D version of Ex6 (Ithmech_2d) with the 3D temperature element (see Ithmech_3d). Afterwards, I have applied the fld3d1tm element and it seems to give same results for steady state and transient case like with the 2D and 3D thermal element.