General formula for dependent probability distributions

Recently I encountered the following problem:

What is the mean distance between two random points on a unit square?

I understand pen and paper methods for solving this exist however I'm interested in a general formula for problems of this type. What I arrived at is the following integral

$$\int_0^1{\int_0^1{\int_0^1{\int_0^1{ \sqrt{(a-b)^2 + (c-d)^2} \space{}da \space{} db \space{}dc \space{}dd}}}}$$

Solving it would be ugly, easily possible with a computer (I get around $0.52$ which if memory serves me matches the non-brute force way). That got me thinking what would be the distribution of the distance between two points on a unit square?

Note: I already calculated it using a simulation and this is the graph, but is there a non-simulation way?:

Ok clearly too complicated to visualise so lets chose a simpler problem for discussion purposes. Lets say the distance between two points on a unit segment (this is just a example, ideally I'm curious for the formula for the probability density of any function of multiple non-uniform probability random variables)

My initial thoughts for this problem is to do some sort line integral under the joint probability density function (which is constant in this case) following the path of the solution of $z=abs(x-y)$ for a given value of $z$ (but to get the formula for all values of $z$. We would be finding the horizontal length of the red surfaces:

We would of course ignore all areas outside the square $0<=x<=1,0<=y<=1$ and do so for all heights of the blue surface. Doing line integrals is analogous to simply reading of the value in a 1-d probability density distribution so we would need a extra formalism regarding the thickness of this curtain under the joint distribution.

I'd expect the formula to be linear which is confirmed by running simulations however I still wonder is this reasoning in any way correct?

After searching for the formula online (which was hard as I have no idea how to specify it) I came across this article on wikipedia:

f(x1, …, xn) shall denote the probability density function of the variables that y depends on, and the dependence shall be y = g(x1, …, xn). Then, the resulting density function is:

$$\int\limits_{y = g(x_1, \cdots, x_n)} \frac{f(x_1,\cdots, x_n)}{\sqrt{\sum_{j=1}^n \frac{\partial g}{\partial x_j}(x_1, \cdots, x_n)^2}} \; dV$$ where the integral is over the entire (n-1)-dimensional solution of the subscripted equation and the symbolic dV must be replaced by a parametrization of this solution for a particular calculation; the variables x1, …, xn are then of course functions of this parametrization.

I get from where the $f()$ came at the top, and I almost get the idea of the parametrization of $dV$ if I try to fit it towards what I believe should be happening above.

What I really don't get is why is the integrand divided by the square root of the laplacian of $g$ (in the line example it was constant so it didn't have a effect on the shape).

Could somebody offer me an explanation for this formula, and most importantly where it comes from?

A recommendation for a book or lectures where to learn more would be the icing on the cake

Sorry for this monster of a question but its been bugging me for a couple days