Note that the maximum initial damage cannot be larger than $D_{max}$. The number of defects $N$ per unit volume of matter
in the range $D_0$ to $D_{max}$ can be calculated by integrating $f(D)$ from $D_0$ to $D_{max}$:

$N = \displaystyle{\int_{D_0}^{D_{max}}} f(D) \mathrm{d} D$

Based on the assumed damage distribution $f(D)$ one can show that
the probability $p$ of having at least one initial defect larger than or equal to $D_0$ in a volume $v$ is:

$p = 1 - \mathrm{e}^{-N \cdot v}$

This probability expression can be used to assign an initial damage level to each integration
point in the model. The damage level is obtained by solving the expression for $D_0$ (given a random number $p$
and an integration point volume $v$).