So, I need to find the center of mass from that trapezoid using double integration knowing the mass function: $\delta(x, y) = 1 + 2x + y$
What I've found so far is the boundaries are $y = -x + 3$, $y = 1$, $x$-axis, $y$-axis
But then, I confused about what is the next step, since I need to use the mass function and use the double integration to find the answer.
How do I find the centroid using double integration and the given mass function?

There is one mistake in the other answer, because the total mass should be computed using the function $\delta$.
Let $M$ be the total mass of the given shape. Then, $$ M = \int_0^1 \left( \int_0^{3-y} \delta(x,y)\, dx\right)dy. $$ First, we shall compute the inner integral: $$ \begin{align*} \int_0^{3-y} \delta(x,y)\, dx &= \int_0^{3-y} (1+2x+y)\, dx \\ &= \left. \frac{(1+2x+y)^2}{4} \right|_0^{3-y}\\ &= \frac{(7-y)^2}{4} - \frac{(1+y)^2}{4} \\ &= \frac{48 - 16y}{4} \\ &= 12 - 4y. \end{align*} $$ So, $$ \begin{align*} M &= \int_0^1 (12-4y)\, dy \\ &= \left. \frac{(12-4y)^2}{-8} \right|_0^1 \\ &= -8 - (-18) \\ &= 10. \end{align*} $$ Hence, the total mass is $M = 10$.
I have written out the rest of the computation as well, so that you can check your calculations.
Next, we want to find the coordinates of the centroid, $(\bar{x},\bar{y})$. To find $\bar{x}$, consider an infinitesimal mass element at the point $(x,y)$. It has mass $\delta(x,y)\, dx\, dy$. We weight this infinitesimal mass by the position coordinate $x$, and then integrate over the whole shape. This value divided by the total mass will give $\bar{x}$. Similarly, for $\bar{y}$, we weight the infinitesimal mass element $\delta(x,y)\, dx\, dy$ by $y$ and perform the integration as before. Mathematically, $$ \begin{align*} M \bar{x} &= \int_0^1 \left( \int_0^{3-y} x (1+2x+y)\, dx \right) dy \\ M \bar{y} &= \int_0^1 \left( \int_0^{3-y} y (1+2x+y)\, dx \right) dy. \end{align*} $$ For $\bar{x}$, the inner integral evaluates as: $$ \begin{align*} \int_0^{3-y} x (1+2x+y)\, dx &= \left. \frac{(1+y)x^2}{2} + \frac{2x^3}{3} \right|_0^{3-y} \\ &= \frac{(1+y)(3-y)^2}{2} + \frac{2(3-y)^3}{3} \\ &= \frac{(3-y)^2(15-y)}{6}. \end{align*} $$ For $\bar{y}$ the inner integral was already computed earlier, and it evaluates to: $$ \int_0^{3-y} y (1+2x+y)\, dx = y(12-4y). $$ So, $$ \begin{align*} M\bar{x} &= \int_0^1 \frac{(3-y)^2(15-y)}{6} \, dy \\ &= \int_0^1 \frac{135 - 99y + 21y^2 - y^3}{6}\, dy \\ &= \left. \frac{45y}{2} - \frac{33y^2}{4} + \frac{7y^3}{6} - \frac{y^4}{24} \right|_0^1 \\ &= \frac{45}{2} - \frac{33}{4} + \frac{7}{6} - \frac{1}{24} \\ &= \frac{123}{8} \end{align*} $$ and $$ \begin{align*} M\bar{y} &= \int_0^1 y(12-4y) \, dy \\ &= \int_0^1 (12y - 4y^2)\, dy \\ &= \left. 6y^2 - \frac{4y^3}{3} \right|_0^1 \\ &= 6 - \frac{4}{3} \\ &= \frac{14}{3}. \end{align*} $$ Therefore, $$ \begin{align*} \bar{x} &= \frac{123}{80} \\ \bar{y} &= \frac{7}{15}. \end{align*} $$