The primary objectives of this paper are to present the construction of bivariate fractal interpolation functions over triangular interpolating domain using the concept of vertex coloring and to propose a double integration formula for the constructed interpolation functions. Unlike the conventional constructions, each vertex in the partition of the triangular region has been assigned a color such that the chromatic number of the partition is 3. A new method for the partitioning of the triangle is proposed with a result concerning the chromatic number of its graph. Following the construction, a formula determining the vertical scaling factor is provided such that the actual double integral coincides with the integral value calculated using fractal theory. Convergence of the proposed method to the actual integral value is proven with sufficient lemmas and theorems. Adequate examples are presented to illustrate the method of construction and to verify the value of double integration.