Let $x\in(0,1)$ be a random variable that follows a truncated normal distribution with density $$ f(x)= \begin{cases} \frac{1}{\sqrt{2\pi\sigma^2}D}\exp\left(-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right) & \text{, if}\:\:0<x<1 \\ 0 & \text{, otherwise} \end{cases} $$ where $D=\Phi\left(\frac{1-\mu}{\sigma}\right)-\Phi\left(-\frac{\mu}{\sigma}\right)$ (let $\Phi$ denote the cdf of the standard normal distribution).
I am interested in evaluating the integral $$ I = \int_{0}^{1}g(x)f(x)\mathrm{d}x, $$ where $g(x)=\log(x)$, $x>0$, is the neperian logarithm, i.e. the integral $$ I = \frac{1}{\sqrt{2\pi\sigma^2}D}\int_{0}^{1}\log(x)\exp\left(-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right)\mathrm{d}x, $$ which, using the substitution $u=\frac{x-\mu}{\sigma}$, is rewritten as follows $$ I = \frac{1}{\sqrt{2\pi}D}\int_{-\frac{\mu}{\sigma}}^{\frac{1-\mu}{\sigma}}\log\left(\sigma u+\mu\right)\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u. $$
I know that this couldn't be given in a closed form. If this makes any difference, I want to use this for a practical application, so if there could be a numerical approximation with some given accuracy, that would suit me, but I would like to avoid it.
First Approach: EDIT: I'm not too fond of this one, the second one seems better to me.I'll write $a = -\frac{\mu}{\sigma}$ and $b = \frac{1-\mu}{\sigma}$ and assume $\mu \neq 0$.
$$I = \frac{1}{\sqrt{2\pi}D}\int_{a}^{b}\log\left(\sigma u+\mu\right)\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u.$$
If $\mu \neq 0$, we can rewrite $\log(\sigma u + \mu) = \log(\frac{1}{\mu}) + \log(1+\frac{\sigma}{\mu}u)$. I'll denote $c := \frac{\sigma}{\mu}$. We can split this up into two integrals $I = I_1 + I_2$, where
$$ I_1 = \frac{1}{\sqrt{2\pi}D} \int_{a}^{b} -\log(\mu) \exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u = \frac{-\log(\mu)}{D}\text{erf}(u)\Big|_{a}^{b},$$ and $$I_2 = \frac{1}{\sqrt{2\pi}D}\int_{a}^{b}\log\left(1+ cu\right)\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u.$$ Here I see two possiblites to continue: Plugging in the taylor expansion for $\log(1+cu)$ or the taylor expansion for $\exp(u)$. I'll begin with the former. Exchanging summation and integral (which i'm not yet sure is legal), and integrating by parts we obtain \begin{align} I_2 &= \frac{1}{\sqrt{2\pi}D}\int_{a}^{b}\sum_{n=1}^{\infty} \frac{(-1)^n}{n} (cu)^n\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u \\ &=\frac{1}{\sqrt{2\pi}D}\sum_{n=1}^{\infty}\int_{a}^{b} \frac{(-c)^n}{n} u^n\exp\left(-\frac{1}{2}u^2\right)\mathrm{d}u\\ &=\frac{1}{D}\sum_{n=1}^{\infty} (-c)^n \left[ \frac{u^n}{n}\text{erf}(u)\Big|_{a}^{b} - \int_{a}^{b}u^{n-1}\text{erf}(u) \mathrm{d}u\right]\\ \end{align}
Looking up the last integral in this integral table (page 5, number 7) , we obtain
\begin{align} I_2 = \frac{1}{D}\sum_{n=1}^{\infty} \frac{(-c)^n}{\sqrt{\pi}(n+1)} \left[ e^{-u^2}\sum_{k=0}^{l-1}\frac{\Gamma(\frac{n}{2}+1)}{\Gamma(\frac{n}{2}-k+1)} u^{n-2k} - (1-j)\Gamma\left(l+\frac{1}{2}\right)\text{erf}(u)\right]_{a}^{b} \end{align}
where $j = 0$ or $j=1$ such that $2l-j = n+1$. Now for the second approach, that yields a much nicer-looking sum.
Second Approach: Starting with
$$I = \frac{1}{\sqrt{2\pi\sigma^2}D}\int_{0}^{1}\log(x)\exp\left(-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right)\mathrm{d}x$$
I'll expand the exponential function: \begin{align} \sqrt{2\pi\sigma}D\cdot I &= \int_{0}^{1} \log(x) \sum_{n=0}^{\infty} \frac{(-1)^n}{2^n \cdot n!}\left(\frac{x-\mu}{\sigma}\right)^{2n} \mathrm{d}x \\ &= \sum_{n=0}^{\infty} \frac{(-1)^n}{\sigma^{2n}\cdot 2^n \cdot n!} \int_{0}^{1} \left(x-\mu\right)^{2n} \log(x) \mathrm{d}x \\ &= \sum_{n=0}^{\infty} \frac{(-1)^n}{\sigma^{2n}\cdot 2^n \cdot n!} \sum_{k=0}^{2n} \binom{2n}{k}(-\mu)^{2n-k} \int_{0}^{1}x^{k} \log(x) \mathrm{d}x \\ \end{align}
Applying the formula $$\int x^{k}\ln x\,dx=\frac{x^{k+1}((k+1)\ln x-1)}{(k+1)^2}$$ we obtain
$$ I = \frac{1}{\sqrt{2\pi\sigma^2}D} \sum_{n=0}^{\infty} \frac{(-1)^{n+1}}{\sigma^{2n}\cdot 2^n \cdot n!} \sum_{k=0}^{2n} \binom{2n}{k}(-\mu)^{2n-k} \frac{1}{(k+1)^2} $$
It's not too hard to see that interchanging integral and sum is allowed using Fubini. We can find yet another representation for $I$: Using the binomial theorem, we see that
$$ \int_{0}^{1}\int_{0}^{1} (xy - \mu)^n \mathrm{d}x\mathrm{d}y = \sum_{k=0}^{2n} \binom{2n}{k}(-\mu)^{2n-k} \frac{1}{(k+1)^2}$$
Thus, we can write \begin{align} I &= \frac{1}{\sqrt{2\pi\sigma^2}D} \sum_{n=0}^{\infty} \frac{(-1)^{n+1}}{\sigma^{2n}\cdot 2^n \cdot n!} \int_{0}^{1}\int_{0}^{1} (xy - \mu)^n \mathrm{d}x\mathrm{d}y \\ &= -\int_{0}^{1} \int_{0}^{1} \frac{1}{\sqrt{2\pi\sigma^2}D} \exp\left(-\frac{1}{2}\left(\frac{xy - \mu}{\sigma}\right)^2\right) \mathrm{d}x\mathrm{d}y \end{align} Now we have a smooth integrand on our domain. Therefore, it should be possible to obtain good numerical results by applying high-degree 2-dimensional tensor product formulas of univariate Gauss-Legendre Quadrature formulas.