The Fokker-Planck partial differential equation (PDE) defines the time evolution of the probability density function of a stochastic process. When the stochastic process is distribution-dependent, the resulting Fokker-Planck PDE is nonlinear. The current work presents a numerical solution to a two-dimensional nonlinear Fokker-Planck PDE using a variational approach. B-Splines with non-uniform knot vectors and finite elements of different orders are used for space discretization, and their convergence is compared. Two types of boundary conditions: absorbing and reflecting are considered in the PDE. The obtained nonlinear system of ordinary differential equations is solved in the time domain by implicit Runge-Kuta methods. The initial condition, which is the Dirac delta function, requires a fine mesh of elements (or knots for the B-Spline discretization) around the non-zero domain of the Dirac delta function. Furthermore, the initial condition requires very small time steps for the time discretization method. Since the Fokker-Planck equation has a diffusion term, the probability density function evolves with time, and the fine mesh of elements (knots) is not required, as it is at the initial time steps. Adaptive mesh generation and adaptive time steps for the space and for the subsequent time discretization of the PDE are used for the numerical solution. Convergence with the order of shape functions is analyzed and presented. The computational time for different adaptive mesh coarsening and adoptive time steps is compared. The numerical experiments are presented using the Heston stochastic local volatility process, which is used to price barrier options. Validation by Monte Carlo simulations is also included.