Do triangular/area coordinates offer computational advantages for numerically integrated triangles?

Isoparametric triangular finite elements can be formulated using "triangular coordinates" $(\zeta_1, \zeta_2, \zeta_3)$, also sometimes called area/areal/barycentric coords, or using "Cartesian coordinates" $(r,s)$ defined on a master element. These two systems can be related to each other without difficulty, and seem to be functionally equivalent.

It seems that several texts and journal papers use triangular coords $(\zeta_i)$ exclusively, particularly those that are more recent and/or mathematically oriented. On the other hand, I have older textbooks and use certain commercial FEA packages that consistently present element formulations with $(r,s)$ coordinates. I have a feeling that this might be due to historical reasons, but perhaps there is more to it.

Provided that the element geometry is restricted (straight sides, and any side nodes are uniformly spaced), I understand that triangular coords may be useful because they lend themselves to analytical integration, and thus avoid potentially costly numerical integration procedures. However, when these conditions are not met, numerical integration must be used. In that case, is there any benefit from using triangular coordinates, conceptual or computational?

Triangular coordinate systems (or tetragonal systems, in 3D) have the advantage of giving all the vertices of the element equal status. That becomes clear in numerical integration procedures, when you look at the positions and weights of the integration points.

The three coordinates are not independent - they should sum to 1 at every point in the element. But usually, preserving the symmetry in the formulas is more useful than losing the independence of the coordinates.

Early papers often picked one vertex (arbitrarily) as "special", and mapped that vertex to the origin of a local system and the adjacent sides along the X and Y axes. For an arbitrary shaped triangular element, this will give different results depending which vertex you put at the origin, which can be confusing when trying to understand why two identical-looking models don't give the same answers!

Before the era of automatic mesh generation and refinement, triangular meshes often had a systematic bias anyway (e.g. start with a rectangular mesh and draw all the diagonals pointing in the same direction) so treating one vertex of the triangle as "special" might have been a good thing to do, in conjunction with that type of mesh. But that's all history now.

Of course you could re-map a "symmetrical" numerical integration procedure onto an "old-school" local coordinate system, if you really wanted to - but it's easier not to bother doing that.

Hardly anyone uses analytic integration for "real world" problems any more - even if the geometry is suitable for it, properties like Youngs Modulus often vary (e.g. with temperature) within a single element. In any case, you can create numerical integration rules that give exact results for the special cases that can be integrated analytically, so why bother to use two different integration methods when one method will handle both situations?

• Re: "For an arbitrary shaped triangular element, this will give different results depending which vertex you put at the origin". For example, are you saying that if we have the 2D (r,s) natural coord sys, with linear shape fcns (N1,N2,N3)=((1-r-s),(r),(s)), and assuming (node1,node2,node3)=((0,0),(1,0),(0,1)), an analysis would have diff't results for any other permutation of the node/shape fcn sequence? Sep 8, 2017 at 23:44
• Also, I would guess that one of the benefits of using analytical integration over the element would be the speedup as compared to quadrature. It may be pointless compared to other factors, but I'm curious: have you seen formulations that proposed to switch from analytical to numerical integration if there is any deviation from the straight-sided restriction? I suppose the goal being that "faster" calcs can be used whenever possible... Sep 8, 2017 at 23:56
• To be honest, I can't think of any current commercial FE software that uses analytic integration. Shaving a few nanoseconds off the execution time is irrelevant compared with a consistent code base, easier quality assurance testing, etc. If you have a simple element (e.g. a linear constant strain triangle) you only need one numerical integration point anyway, so the computational cost of the integration is mainly computing the Jacobian of transformation from real-world coordinates to element coordinates (in other words, what is the real-world area of the element), which is needed either way. Sep 9, 2017 at 0:34
• Re "different results", it depends on where the integration points are in the element coordinate system. If you do analytical integration correctly, it would make no difference. But given a 90-45-45 triangle in the element coordinate system some the "obvious" numerical integration schemes (at least, back in the 1960s!) were not symmetrical with respect to the vertices of the triangle. Sep 9, 2017 at 0:38
• For large models the cost is dominated by the equation solution, not the numerical integration. The economics of this has changed a lot, when "large" now means maybe 1,000,000 degrees of freedom, compared with maybe 10,000 DOF 50 years ago, and computing speeds have increased by more than 1,000 times over that half-century - not to mention new equation solution algorithms, etc. Sep 9, 2017 at 0:42