Given the problem $$-\nabla u + u = f$$ $$ n\cdot\nabla u = g \quad\text{on} \quad \Gamma$$ I can show the discretization given through the Galerkin formulation is $Au=b$ where $$ A = \int_\Omega \nabla\varphi_i\cdot\nabla\varphi_jdx + \int_\Omega \varphi_i\varphi_jdx$$ $$ b= \int_\Omega f \varphi_idx + \int_\Gamma g\varphi_ids $$ and the basis functions are linear with form $\varphi_i = c_1+c_2x+c_3y$, satisfying the condition that $$\varphi_i(x_j) = \delta_{ij}$$ For $A$, I can find the first term's elemental stiffness matrix through $$A = \int_\Omega (c_{2,i}c_{2,j} + c_{3,i}c_{3,j})dx = Area(c_{2,i}c_{2,j} + c_{3,i}c_{3,j})$$
I'm trying to calculate the elemental stiffness matrix for $$\int_\Omega\varphi_i\varphi_jdx$$ but getting stuck as I'm not sure what do with $c_3$ terms
Since $\varphi_i \varphi_j$ is a polynomial of degree 2, you can just use a Gauss quadrature formula (if the formula you choose has degree $\ge 2$, it will yield the exact value for the integral). You can either compute the quadrature nodes form each element, or compute all the integrals in a reference element, using a change of variables.