Yeah, it works fine using the penalty method, although for me personally lagm is a far more elegant method and shows better results in serial solution.

I would like to comment on a little offtopic bug:
Feap offers the possibility to include a rigid contact surface into the calculation.
On p. 173 of the FEAP 8.6 User manual it is written, that the surface normal for a cartesian plane is either specified by positive or negative coordinate direction.
However, while using negative direction say -3 for negative z-direction, Feap prompts an error: "**ERROR** No user rigid surface -3" and stops.
I was able to circumvent this error by commentic the plstop in line 77 in $FEAPHOME8_6/contact/ntrnd/cruser.f and achieved reasonable results.
For me it is quite hard to find the origin of the error.
Maybe you could look at the bug and fix it for a future release.
P.S.: A demo input file is attached.