Are there any two-dimensional quadrature that only uses the values at the vertices of triangles?

Solution 1:

One-dimensional

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:

• Understanding Galerkin method of weighted residuals

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:

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.

Two-dimensional

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:

• What is the difference between Finite Difference Methods, Finite Element Methods and Finite Volume Methods for solving PDEs?

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:

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 .

Three-dimensional

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:

• Presentation at Woudschoten 1996
• On solving a Cosmic Ray equation
• Resistor Models for Diffusion in 3-D
• Recipe for Convection & Diffusion in 3-D
• CR pressure with AVS