After reading on my notes about two dimensional quadrature rules, I noticed that there weren't any Gaussian quadrature rules (or any in general) that only used values at the vertices of the triangles. My question arises when thinking about a triangular mesh wherein the values of a function are only known at the mesh points, and you want to approximate the values of the function over the domain.

What would be the best way to go about doing something like this? Interpolating from the values at the nodes to the quadrature points (like the centroid) or is there a way to derive a quadrature rule that uses only the vertices?

  • 1
    $\begingroup$ Generally, quadrature methods try to avoid the boundary of the domain, so that they don't evaluate the function at singularities at the endpoints. You'll probably have to write down your own quadrature method to do what you want. $\endgroup$ – user14717 May 3 '18 at 1:27
  • $\begingroup$ I see, @user14717 can you point me somewhere where the computation of 2D quadratures on triangles are discussed? $\endgroup$ – Gonate May 4 '18 at 17:12
  • 1
    $\begingroup$ Because my expertise has become a bit rusty from time tot time, I've made some wrong statements in the first version of my answer below. Which have been corrected now. But I'm not finished yet: please be patient. $\endgroup$ – Han de Bruijn May 8 '18 at 16:42
  • 1
    $\begingroup$ I've completed my answer with 2-D and 3-D. Hope you appreciate. $\endgroup$ – Han de Bruijn May 8 '18 at 20:29


Let's go back to the basics. Compute the integral at the interval $[-1/2,+1/2]$ of a function $f(\xi)$ numerically, by employing a 2-point Gauss-Legendre integration: $$ \int_{-1/2}^{+1/2} f(\xi)\,d\xi = w_1.f(\xi_1)+w_2.f(\xi_2) $$ Determine weights $(w_1,w_2)$ and locations $(\xi_1,\xi_2)$ of the so-called integration points by the requirement that a polynomial of degree $3$ is integrated exactly. Four equations can be set up, as follows: $$ 1 = \int_{-1/2}^{+1/2} 1\,d\xi = w_1.1 + w_2.1 \quad \mbox{(1)} \\ 0 = \int_{-1/2}^{+1/2} \xi\,d\xi = w_1.\xi_1 + w_2.\xi_2 \quad \mbox{(2)}\\ \frac{1}{12} = \int_{-1/2}^{+1/2} \xi^2\,d\xi = w_1.\xi_1^2 + w_2.\xi_2^2 \quad \mbox{(3)}\\ 0 = \int_{-1/2}^{+1/2} \xi^3\,d\xi = w_1.\xi_1^3 + w_2.\xi_2^3 \quad \mbox{(4)} $$ These four independent equations are sufficient to determine two weights $(w_1,w_2)$ and two locations $(\xi_1,\xi_2)$ . Now the two equations (2),(4) are easily satisfied by adopting the following equalities: $w_1 = w_2$ , $\xi_1 = - \xi_2$ . Substitution into equation (1) gives : $w_1 = w_2 = 1/2$ . Substitute this into the remaining equation (3), giving: $1/12 = \xi_2^2$ . Hence $\xi_2 = \sqrt{3}/6$ . Summarizing: $$ (w_1,w_2) = (\frac{1}{2},\frac{1}{2}) \quad \mbox{and} \quad (\xi_1,\xi_2) = (-\frac{\sqrt{3}}{6},+\frac{\sqrt{3}}{6}) $$ In Finite Element culture, sometimes reduced integration is employed, meaning in our case that $f(\xi)$ shall be integrated as if it were only a polynomial of degree $1$, which is linear. In this case we have only the equations (1) and (2) that must be satisfied: $1 = w_1 + w_2$ and $0 = w_1.\xi_1 + w_2.\xi_2$ . There are a myriad solutions for this to accomplish. Let's assume for reasons of symmetry that $w_1 = w_2 = 1/2$ . Then two extreme cases are $\xi_1=\xi_2=0$ and $(\xi_1,\xi_2) = (-1/2,+1/2)$ .
The first case is also covered by a one-point Gaussian quadrature, which is the better way to do it: $$ 1 = \int_{-1/2}^{+1/2} 1\,d\xi = w_1 \quad ; \quad 0 = \int_{-1/2}^{+1/2} \xi\,d\xi = w_1\xi_1 \quad \Longrightarrow \quad (w_1,\xi_1) = (1,0) $$ The second case is integration at the nodal points, according to the OP's wishful thinking.

Let's work out the above for a specific problem that has been solved recently at MSE:

In this answer we find the following expressions for the finite element matrix (look it up): $$ E_{0,0}^{(i)} = E_{1,1}^{(i)} = 1/(x_{i+1}-x_i)+(x_{i+1}-x_i)\,p^2/3 \\ E_{0,1}^{(i)} = E_{1,0}^{(i)} = -1/(x_{i+1}-x_i)+(x_{i+1}-x_i)\,p^2/6 $$ If we hadn't worked out the integrals, then this would have been: $$ E_{0,0}^{(i)} = E_{1,1}^{(i)} = 1/(x_{i+1}-x_i)+(x_{i+1}-x_i)\,p^2 \int_{-1/2}^{+1/2} \left(\frac{1}{2}\pm\xi\right)^2 d\xi \\ E_{0,1}^{(i)} = E_{1,0}^{(i)} = -1/(x_{i+1}-x_i)+(x_{i+1}-x_i)\,p^2 \int_{-1/2}^{+1/2} \left(\frac{1}{4}-\xi^2\right) d\xi $$ The integrals can be worked out with integration points as well. First take the common Gauss points: $$ \int_{-1/2}^{+1/2} \left(\frac{1}{2}\pm\xi\right)^2 d\xi = \frac{1}{2}\left(\frac{1}{2}-\sqrt{3}/6\right)^2+\frac{1}{2}\left(\frac{1}{2}+\sqrt{3}/6\right)^2 = \frac{1}{4}+\frac{1}{12} = \frac{1}{3} \\ \int_{-1/2}^{+1/2} \left(\frac{1}{4}-\xi^2\right) d\xi = 2\times \frac{1}{2}\left[\frac{1}{4}-\left(\pm\sqrt{3}/6\right)^2\right] = \frac{1}{4} - \frac{1}{12} = \frac{1}{6} $$ Exactly as it should, of course. Now we do reduced integration for the two extreme cases. First case is center of element $\xi = 0$ : $$ \int_{-1/2}^{+1/2} \left(\frac{1}{2}\pm \xi\right)^2 d\xi = 1.\left(\frac{1}{2}\pm 0\right)^2 = \frac{1}{4} \\ \int_{-1/2}^{+1/2} \left(\frac{1}{4}-\xi^2\right) d\xi = \frac{1}{4} - 0^2 = \frac{1}{4} $$ Last but not least the case that is wanted by the OP, reduced integration at the nodal points of the element $(\xi_1,\xi_2) = (-1/2,+1/2)$ : $$ \int_{-1/2}^{+1/2} \left(\frac{1}{2}\pm \xi\right)^2 d\xi = \frac{1}{2}\left(\frac{1}{2}-\frac{1}{2}\right)^2 + \frac{1}{2}\left(\frac{1}{2}+\frac{1}{2}\right)^2 = \frac{1}{2} \\ \int_{-1/2}^{+1/2} \left(\frac{1}{4}-\xi^2\right) d\xi = 2\times \frac{1}{2}\left[\frac{1}{4}-\left(\pm\frac{1}{2}\right)^2\right] = 0 $$ Thus we now have for the finite element matrix at hand: $$ E_{0,0}^{(i)} = E_{1,1}^{(i)} = 1/(x_{i+1}-x_i)+(x_{i+1}-x_i)\,p^2 \times \begin{cases} 1/3 \\ 1/4 \\ 1/2 \end{cases} \\ E_{0,1}^{(i)} = E_{1,0}^{(i)} = -1/(x_{i+1}-x_i)+(x_{i+1}-x_i)\,p^2 \times \begin{cases} 1/6 \\ 1/4 \\ 0 \end{cases} $$ Question is: which one is the best? All these integration schemes seem to be equally good at first sight. At first sight, because there is another issue to be covered : stability / robustness .

The following book is downloadable from the internet. It's about Finite Volume methods, but it contains important guidelines for the Finite Element method as well.

  • Suhas V. Patankar, Numerical Heat Transfer and Fluid Flow.
According to one of Patankar's "Four Basic Rules" (at page 37), "the rule of positive coefficients", off-diagonal terms $E_{i,j}$ with $i \ne j$ must be less than zero. With our 1-D problem this means a restriction on the mesh size for common and one-point integration (same colors are employed in the picture below): $$ E_{0,1}^{(i)} = E_{1,0}^{(i)} = -1/(x_{i+1}-x_i)+(x_{i+1}-x_i)\,p^2 \times C < 0 \\ \Longleftrightarrow \quad C\,\left[\,p\,(x_{i+1}-x_i)\,\right]^2 < 1 \quad \mbox{with} \quad C = \begin{cases} \color{blue}{1/6} \\ \color{red}{1/4} \\ \color{green}{0} \end{cases} $$ The above condition becomes increasingly difficult to fulfill when $\,p\,$ is large. The only exception being $\,C=\color{green}{0}$ , which happens to be : integration (points) at the element vertices !
This is seen in the following picture for our constant $p=100$ and number of mesh points $N=20$. Mind the different colors:
enter image description here
Thus we may conclude that the case with $\,\color{red}{C = 1/4}\,$ (i.e. centroid integration) is worst and the case with $\color{green}{C = 0}$ (i.e vertex integration) is best.
Especially in higher dimensional problems, the coefficients of equations may vary wildly from place to place. So robustness (i.e. unconditional stability) is not a luxury, certainly not in 2-D and 3-D.


But there is more. Alternative integration schemes in FEM are sort of a stargate for incorporating knowledge from other numerical methods, such as the Finite Volume Method. As is exemplified here for example:

In that answer it has not been mentioned that integration points at a quadrilateral finite element are actually taken at the nodal points, i.e. the vertices, as has been suggested in the OP's question:
enter image description here
Resulting in a splitting of the quadrilateral in four linear triangles:

Finally, these triangles comprise an equivalent finite volume method for the diffusion problem in a finite element context. A much more complete account of this theory is found elsewhere, being too large to fit into the margins of MSE : 2-D Elementary Substructures .


The higher the dimension, the more difficult it is to formulate a concise answer. But it is the same story all over the place: integration points at the vertices (of a brick in 3-D) are the best option to guarantee unconditional stability of finite element schemes. Where it should be emphasized that numerical schemes don't have to be the same for different terms of the partial differential equation. The result of the research below has been the publication of a paper by Horst Fichtner, this author and others. The title of the paper is Longitudinal gradients of the distribution of anomalous cosmic rays in the outer heliosphere. Some (numerical and graphical) details are found in the following references:

  • $\begingroup$ Thanks, I've been looking through all the references included in your answer these past days. The first part was all I really needed, the problem wasn't that difficult, but I plan to read through all the extra material included. $\endgroup$ – Gonate May 11 '18 at 13:40
  • 1
    $\begingroup$ @Gonate: Maybe another nice reference : Circuits on an arbitrary triangular lattice . $\endgroup$ – Han de Bruijn May 11 '18 at 20:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.