We evaluate the integral expressions for the matrices (12), (14), and (16) numerically using Gaussian quadrature . We shall explain the computation of the element stiffness matrix; the computation of the mass and damping matrices follow suit. Assuming the element's parametric domain is , the expression for entry of the stiffness matrix of a D-NURBS surface element takes the integral form
where, according to (16),
Here, the are the columns of the Jacobian matrix for the D-NURBS surface element.
We apply Gaussian quadrature to compute the above integral approximately. The integral is obtained by applying Gaussian quadrature on the 1-D interval twice. Given integer and , we can find Gauss weights , and abscissas , in two directions of the parametric domain such that can be approximated by ()
We apply the de Boor algorithm  to evaluate .
Generally speaking, for integrands that are polynomial of degree 2N-1 or less, Gaussian quadrature evaluates the integral exactly with N weights and abscissas. For D-NURBS, is not polynomial unless the model is reduced to a B-spline. In our system, we choose and to be integers between 4 and 7. Our experiments reveal that matrices computed in this way lead to stable, convergent solutions.