Hi,
I ve written a program for solving Reynolds equation using FEA.
I can give you guidelines :
1-/ mesh the 2D domain of sliding (I used triangle 3 nodes)
and apply boudary conditions (p=0 all around the domain for starting)
2-/ build the element matrix, (ref: Huebner "FEM for engineers", Reddi, and so on). For the film thickness, use numerical integration for h and h^3, very simple for linear shape function
3-/ assemble the general matrix which is sym. and positive definite, that is no problem for solving.
4-/ Post process the pressures to obtain velocity, shear stress, load and so on.
Next you can add enhancements such as error estimation, embedded sparse solver (use some PCG), refine polynomial for the recovery of second order derivatives.
I remains avalaible to help you in that subject.
Regards