I assume you're trying to solve an equation that looks like:
\begin{align}
-\nabla \cdot (a(x)\nabla{u}) = f,
\end{align}
for $x$ in some domain $\Omega$, although the same approach would be fine (for a residual evaluation, anyway) if $a$ were also a function of $u$.
The stiffness matrix will take the form
\begin{align}
A_{ij} = \int_{\Omega}a(x)\nabla\varphi_{j}(x)\cdot\nabla\varphi_{i}(x)\,\mathrm{d}x.
\end{align}
For $P_{1}$ triangular elements, $\nabla\varphi_{j}$ will be constant over each element $K$, for each $j$. Letting $\mathcal{K}$ be the mesh, the elements of the stiffness matrix become
\begin{align}
A_{ij} = \sum_{K \in \mathcal{K}}\nabla\varphi_{j}^{K}\cdot\nabla\varphi_{i}^{K}\int_{K}a(x)\,\mathrm{d}x,
\end{align}
where all I've done is factor out the constant gradients over each element (these gradients will differ from element to element, which is why they are also indexed over $K$). If $K$ has vertices $N_{1}$, $N_{2}$, and $N_{3}$, then
\begin{align}
\int_{K}a(x)\,\mathrm{d}x = a\left(\frac{N_{1} + N_{2} + N_{3}}{3}\right)|K| + O(h^{3}),
\end{align}
where $h$ is the mesh "size", and $|K|$ is the area of element $K$. This approximation is also called the center of gravity rule.
In practice, you'd take the local stiffness matrix for each element (assuming a diffusivity of one), and multiply it by the variable diffusivity $a$ evaluated at the center of gravity of the element.