$K u$ equals the internal forces only in the linear case. The tangent stiffness matrix, $K$, in a nonlinear problem is normally used in a Newton-Raphson algorithm to calculate updates to the displacement vector as follows:
$$K \Delta u = f - f_{internal}$$ $$ u_{i+1} = u_i + \Delta u$$
The vector of internal forces, $f_{internal}$ must be calculated from the nonlinear element equations. The iteration continues until $u$ converges.
In a nonlinear problem, both the tangent $K$ and the $f_{internal}$ vector are functions of the displacements. One consequence of this is that to prescribe displacements at selected nodes, the displacement vector passed into the functions that calculate these must contain the prescribed values. The nonlinear solution algorithm typically begins with a solution vector of all zeros. Instead, selected entries can be set to the prescribed values so that the internal forces will be calculated correctly. If VegaFEM explicitly eliminates global equations for a constrained degree of freedom, it also should do that for prescribed displacements. Those might be the only changes to the code needed to prescribe non-zero instead of zero displacements.