For transient problem, we can
 

For transient problem, we can achieve the solution as following:
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == -fvc::grad(p));

Then, if I want to solve a arbitrary Poisson's equation or elliptic equation in my solver code, I can do it with the similar code lines shown above?

Thanks.
BTW, I have successfully include VTK libs in my OpenFOAM solver. Thanks.