This paper concerns with analytical integration of trivariate polynomials over linear polyhedra in Euclidean three-dimensional space. The volume integration of trivariate polynomials over linear polyhedra is computed as sum of surface integrals in R3 on application of the well known Gauss's divergence theorem and by using triangulation of the linear polyhedral boundary. The surface integrals in R3 over an arbitrary triangle are connected to surface integrals of bivariate polynomials in R2. The surface integrals in R2 over a simple polygon or over an arbitrary triangle are computed by two different approaches. The first algorithm is obtained by transforming the surface integrals in R2 into a sum of line integrals in a one-parameter space, while the second algorithm is obtained by transforming the surface integrals in R2 over an arbitrary triangle into a parametric double integral over a unit triangle. It is shown that the volume integration of trivariate polynomials over linear polyhedra can be obtained as a sum of surface integrals of bivariate polynomials in R2. The computation of surface integrals is proposed in the beginning of this paper and these are contained in Lemmas 1–6. These algorithms (Lemmas 1–6) and the theorem on volume integration are then followed by an example for which the detailed computational scheme has been explained. The symbolic integration formulas presented in this paper may lead to an easy and systematic incorporation of global properties of solid objects, for example, the volume, centre of mass, moments of inertia etc., required in engineering design processes.