How to approximate the integral of the product between two hat functions? Consider the partition of the time interval $[0,T]$ in $n$ equispaced subintervals, and consider the family $(\phi_i)_{i=0,...,n}$ of hat functions defined on $[0,T]$ (see Appendix A below).

Objective is to approximate $\int_0^t \mu(s,t) X_s \,ds$ so that the
  approximation will have the form $X^Tv$, where $X$ and $v$ are two
  vectors.

The functions $\mu(s,t)$ and $X_s$ can be approximated by (see Appendix B below) 
$$
   X_s \approx X^T\Phi(s) = \Phi^T(s)X
$$
$$
   \mu(s,t) \approx \Phi^T(s) M \Phi(t) = \Phi^T(t) M^T \Phi(s)
$$
where $\Phi(t) = [\phi_0(t),...,\phi_n(t)]^T$, $X$ is a vector of length $n+1$ and $M$ is a square matrix of size $n+1$. Hence
$$
   \int_0^t \mu(s,t) X_s \,ds \approx X^T\Bigg(\int_0^t \Phi(s)\Phi^T(s) ds\Bigg) M\Phi(t)
$$
A property of the hat functions is that $\int_0^t \Phi(s) \,ds \approx P\Phi(t)$, where $P$ is a square matrix of size $n+1$ called operational matrix of integration.
Following the technique used in this article, denote with $M_i$ and $P_i$ the $i$th rows of matrices $M$ and $P$, then
\begin{align}\tag1
    \color{blue}{\Bigg(\int_0^t \Phi(s)\Phi^T(s) ds\Bigg) M\Phi(t)}
    &{}\color{blue}{\approx \begin{pmatrix}P_0\Phi(t)M_0\Phi(t) \\ \vdots \\ P_n\Phi(t)M_n\Phi(t)\end{pmatrix}} \\[0.5em]
    &{}\color{blue}{\approx \begin{pmatrix}P_0 \text{diag}(M_0) \\ \vdots \\ P_n \text{diag}(M_n)\end{pmatrix} \Phi(t)}\\[0.5em]
    &{}\color{blue}{=  P \odot M \Phi(t)}
\end{align}
where $\odot$ is the element-wise product of matrices.
I don't understand how $(1)$ is computed, maybe we have to use the fact that $\Phi(t)\Phi^T(t) \approx \text{diag}(\Phi(t))$ ?

EDIT 1
It seems that $\int_0^t \Phi(s)\Phi^T(s) ds \approx P$, but if this was the case, then how to explain the use of the element-wise product instead of the regular product?

EDIT 2
I'm trying to use the relation $\Phi(t)\Phi^T(t)v = \text{diag}(v)\Phi(t)$, with $v$ a vector. Write $M$ as a column vector whose elements are the rows of $M$, then
$$
   \int_0^t \Phi(s)\Phi^T(s)M ds = \int_0^t \Phi(s)\Phi^T(s) \begin{pmatrix}M_0 \\ \vdots \\ M_n \end{pmatrix} ds \approx \int_0^t \text{diag}(M_0,...,M_n)\Phi(s)ds = \text{diag}(M_0,...,M_n) \int_0^t \Phi(s)ds \approx \text{diag}(M_0,...,M_n) P\Phi(t)
$$
Then
$$
   \color{blue}{\Bigg(\int_0^t \Phi(s)\Phi^T(s) ds\Bigg) M\Phi(t)} \approx \text{diag}(M_0,...,M_n) P\Phi(t)\Phi(t)
$$
but the product $\Phi(t)\Phi(t)$ is defined only if pointwise, moreover I don't think $\text{diag}(M_0,...,M_n)$ is defined.

Appendix A
Consider the partition of the time interval $[0,T]$ in $n$ equispaced subintervals of length $h=T/n$. The family of $n+1$ hat functions on $[0, T]$ is defined as
$$
    \phi_0(t) =
    \begin{cases}
        \frac{h-t}{h},  &0\le t\le h \\
        0,              &\text{ otherwise}
    \end{cases}\quad
    \phi_n(t) =
    \begin{cases}
        \frac{t-(T-h)}{h},  &T-h\le t\le T \\
        0,                  &\text{ otherwise}
    \end{cases}
$$
$$
    \phi_i(t) =
    \begin{cases}
        \frac{t-(i-1)h}{h}, &(i-1)h\le t\le ih \\ \frac{(i+1)h-t}{h}, &ih\le t\le (i+1)h \\
        0,                  &\text{otherwise}
    \end{cases} \quad\implies\quad
    \begin{aligned}
        &\phi_i(jh) =
        \begin{cases}
            1, &i=j \\
            0, &i\ne j
        \end{cases} \\
        &\phi_i(t)\phi_j(t) = 0, \text{if } |i-j|\ge2
    \end{aligned}
$$
with $i=\overline{1,n-1}$.

Appendix B
A function $f$ on $[0,T]$ can be approximated as
$$
    f(t) \approx \sum_{i=0}^n f_i\phi_i(t) = F^T\Phi(t) = \Phi^T(t)F
$$
where $F=[f_0,...,f_n]^T,\ f_i = f(ih),\ i=\overline{0,n}$ and $\Phi(t) = [\phi_0(t),...,\phi_n(t)]^T$.
A function $g$ on $[0,T]\times [0,T]$ can be approximated as
$$
    g(s,t) \approx \Phi^T(s) \Lambda \Phi(t)
$$
where $\Lambda$ is a matrix given by $g_{ij} = g(ih,jh),\ i,j = \overline{0,n}$.

MH Heydari, MR Hooshmandasl, FM Maleek Ghaini, C. Cattani: "A computational method for solving stochastic Itô–Volterra integral equations based on stochastic operational matrix for generalized hat basis functions", J. Comput. Phys. 270, 402-415. doi:10.1016/j.jcp.2014.03.064
 A: Using the results of the article and this post,
\begin{aligned}
A(t) = \int_0^t \Phi(\tau)\Phi(\tau)^T \text d \tau &\simeq \int_0^t \text{diag}\,\Phi(\tau)\, \text d \tau \\
&= \text{diag} \left(\int_0^t \Phi(\tau)\, \text d \tau \right) \\
&\simeq \text{diag}\big( P \Phi(t) \big) \, .
\end{aligned}
Thus, the diagonal entries of the matrix $A$ are the scalar products $P_i\Phi$, where $P_i$ is the $i$th row of $P$. The vector $M\Phi$ may be viewed as the column vector of the scalar products $M_i\Phi$, where $M_i$ is the $i$th row of $M$. Thus, the product $AM\Phi$ is the column vector of scalars $P_i\Phi\, M_i\Phi$ with $i = 0, \dots , n$. Now, let us "expand the entries" of $AM\Phi$ by the hat functions. If we set
$$
[AM\Phi]_i = P_i\Phi\, M_i\Phi \simeq \sum_{j=0}^n a_{ij} \Phi_j \, ,
$$
then evaluation at the grid nodes $t = \ell h$ gives $P_{i\ell} M_{i\ell} = a_{i\ell}$. Thus, the final approximation
$$
[AM\Phi]_i \simeq \sum_{j=0}^n P_{ij} M_{ij} \Phi_j = \sum_{j=0}^n\, [P \odot M]_{ij} \Phi_j
$$
is obtained. QED.
