Standard literature in finite elements hardly addresses the convexity issue. There is no mention of it in the encyclopedic references in the two volumes of Zienkiewicz and Taylor [9]. Wachspress addresses the condition from a projective geometry point of view. He develops concavity with curved sides. His general procedure [10] for n-gon is not applicable since the curve through the external intersection points intersects the element domain. It is then prudent to conclude that there has not been any consistent method available to generate the shape functions of a concave polygon. Recently, Dasgupta utilized the two-dimensional version of the aforementioned code in the form shapeRelation[f,s,{x,y},{x0,y0}] to generate the second, third, and the fourth shape functions from the s which pertained to the convex end.

The following code calls the shape function generator, p is the list of the nodal coordinates, and the physical domain is represented by the x-y frame. An example with the following nodes:

is presented here.

The linearity condition along the straight boundary segments of the calculated shape functions is verified:

In the matrix, rows indicate the shape function and columns indicate corresponding sides of the quadrilateral.

The original shape function is displayed with unit value at the concave end and with zero values at the remaining ends.

The following three shape functions are calculated from the linear interpolating condition: f={1,x,y}.

The second shape function below is associated with the bottom-most node.

The third interpolant shown below corresponds to the node facing the concave end.

The following interpolant is associated with the top-most node. Very low values can be observed in the bottom region of the element.

The following analytical forms of the shape functions can be used to generate the element stiffness matrix using the exact integration Mathematica code presented at IAS in 1998, in Murmansk, Russia.

Observe that the Padé form is automatically captured.

Let us consider a quadrilateral with the following nodes.

Let the shape function associated with {x1,y1} be n1[x,y] in the x-y coordinate frame. The three other shape functions, n2[x,y], n3[x,y], and n4[x,y] were obtained from:

A remarkable result is that if n1[x,y] satisfies the elliptic condition , then the same partial differential equation, for arbitrary , is satisfied by the three computed shape functions as verified below.