This course will continuously refer to . That book is highly recommended as a reference. Detailed handwritten notes of important derivations for stressRefine are also available here.
For conventional finite elements the basis functions used to represent displacement are continuous at interfaces between elements as long as the adjacent elements have the same nodes at the interface. These functions are often called “shape functions” because they’re the same functions used to represent the element mapping from natural coordinates to physical coordinates, which determines the element’s physical shape. Here will make a distinction. “Basis functions” is used to describe functions representing displacement, while “mapping functions” (interchangeable with “shape functions”) are used to describe functions representing the element mapping.
Before I get to how continuity of basis functions, let me describe a slight tweak I made in stressRefine to the classic heirarchical basis functions developed by Szabo and Babusca .
Variation used in stressRefine:
In the classic functions the only nodes are at the element corners and the edge functions start at polynomial order 2. In the stressRefine functions, the nodal functions are the same as the conventional finite element quadratic functions, which requires the introduction of mid-edge nodes. The edge functions start at polynomial order 3. The advantage of this change is that if the elements is a polynomial order 2, its is identical to a conventional quadratic isoparametric element. A possible disadvantage is if this would degrade the element stiffness conditioning. (The element stiffness matrices are assembled into a linear system Ku = f. If K is not well-conditioned, slight errors in f can cause large errors in u). But I was able to prove numerically that the element stiffness matrices are as well conditioned with the modified functions as with the classic ones, at least as high as polynomial order 8.
So for the stressRefine functions, continuity across element interfaces is assured if the adjacent elements have the same nodes, up until polynomial order 2. Additional work is needed to assure continuity at shared edges and faces for polynomial orders higher than 2.
To show why this is so for the case of edges, consider the odd-numbered polynomial functions for polynomial orders 3,5,7, … In the figure shown (for a 2D mesh), the upper element has the shared edge defined between corner 1 and 2 based on the local element numbering, so it, while the lower elements has the shared edge defined between corners 2 and 3, so the edge runs backwards from the other element. This means that the local definition of the p3 (short for polynomial order 3) function is different, and for the lower element that function has to be multiplied by -1 to assure continuity.
The situation is more complicated for faces. Consider an exploded view of two adjacent elements that share a face with global nodes 10, 11, 12:
But for the left element these correspond to local nodes 1,2,3 while on the right element they correspond to local nodes 3,1,2. So for the two elements the natural coordinates on the face rf, sf are defined differently, and the basis functions are misoriented on the right face with respect to the left.
This brings up an interesting historical aside about Mechanica. My colleague and cofounder Christos Katsis handled everything related to coding up the basis functions. He introduced the concept that there is a global face, and local faces can be rotated (or even reflected) compared to it. There is then a coordinate transformation to make sure the local face is compatible with the global face. It is further complicated by the fact that we need to relate elements natural coordinates r,s,t to face natural coordinates rf, sf. So there was a 3×2 transformation relating r,s,t to the local face, and a second 2×2 transformation relating the local face to the global face. None of this is computationally intensive, but it was tricky to code. Christos took care of all of this while we were in full-blown early start-up mode, with tons of distractions. A major feat of concentration!
When I started to develop stressRefine I was not looking forward to the equivalent of all of that, but fortunately found there’s an easier way to handle it. There was a hint in Szabo and Babusca’s textbook , when describing implementation of the hierarchical basis functions. When discussing the 2D basis functions for a triangle, instead of referring to local coordinates r and s they used l2 – l1 and 2l3-1 where l1,l2, and l3 are the area coordinates of the face. From the definition of the area coordinates, r = l2 – l1 and s = 2l3-1. I realized that they had done this because it guarantees continuity of basis functions for elements sharing the same face as long as the nodes of the face are numbered consistently. So we just need to introduce the concept of how the local node numbering of the face relates to the node numbering of the global face.
The global face is defined by the first element that owns that faces. So in the example above, the left hand element is encountered first, and the global node numbering of the face is local node 1 = global node 10, local node 2 = global node 11, and local node 3 = global node 12. Now when we encounter the face again in the right hand element, the local nodes of the face are 3,1,2. When we work with the basis functions of the face, compatibility is assure if we use rf =l1– l3 and sf = 2l3-1 when referring to the face for the second element. This is for 2D, It turns out in generalizes to 3D if we use the same idea with volume coordinates instead of area coordinate. This also takes care of relating element r,s,t to face rf, sf.
I’ll show how this is coded up below. In short, we introduce a variable gno for “global node order”, so for the example face on the right hand element gno(1) = 3, gno(2) = 1, and gno(3) = 2.
All of this works for quadrilaterals in 2D, and bricks and wedges in 3D as well, But for bricks, for example, we define rf, sf using the linear mapping functions of the brick N1 through N8 instead of volume coordinates.
Details of Basis functions for a tet
Nodal functions: These have the value 1 at one node (corner or midedge) and 0 at all others, so they are the same as the quadratic shape functions for a 10 noded tet . for example N1 = (2L1 -1)L1 is a typical function for a corner node, N5 = 4L1L2 is a typical function for a midedge node.
The edge, face, and volume functions are the same as described in  except the edge functions start at p3.
Edge Functions: These are nonzero on one element edge and zero on all others. They start at p3 because the quadratic edge functions are the same as nodal midedge functions. They are the integrals of Legendre polynomials shown here blended with volume coordinates so they are zeroes on all other edges:
bn = 4LILJφp(re)/(1-re2) where φp is the integral of Legendre polynomial over the edge for polynomial number p, re = LJ – LI is the natural coordinate on the edge, and n is the basis function number of the element. The term LILJ with the volume functions assures this reduces to the nonzero 1D function φp on the edge with corners I, J, and is 0 on all other edges, and zero on all faces except those with corner I,J.
Face Functions: These are nonzero interior to one face, zero on all edges, and blend to 0 on all other faces.
bn = 4LILJLkPpr(rf)Pps(sf) where Ppr is the Legendre polynomial in the rf direction and pr is the polynomial number in that direction, Pps is the Legendre polynomial in the sf direction and ps is the polynomial number in that direction, rf = LJ – LI and Sf = 2Lk-1 are the natural coordinates on the face, and n is the basis function number of the element The term LILJLk with the volume functions assures this is nonzero interior to the face with corners I, J,K, and is 0 on all edges,and all other faces,
Volume functions: These are 0 on all edges and faces, and nonzero in the interior to the element.
bn = LILJLkLlPpr(r)Pps(s)Ppt(t) where Ppr is the Legendre polynomial in the r direction and pr is the polynomial number in that direction, and similarly for s,t. The volume functions assure that these are 0 on all edges and faces.
Wedges and Bricks
Similar techniques are used to blend the 1D functions wedges, and bricks, For calculating element stiffnesses, we need the derivatives of the basis functions with respect to the element natural coordinates r,s,t. This is a little messy but not difficult to code once you have the formulae derived. The basis functions and their derivatives for tets, bricks, and wedges are described in my notes here. Note that we do not have to actually evaluate the integrals of Legendre polynomials, there are regression formulae for computing them described in .
Use of “global node order”: I described assuring continuity of the basis functions above. For edges we just need to correct if the local edge runs backwards from the global edge. A direction variable, which is +/- 1 is assigned to each edge. The basis function is multiplied by “direction” to assure continuity (the same thing could be achieved by switching I and J in the edge basis function if the edge runs backwards, but the direction flag is easier. For the face functions of a triangular face we need to use the “global node order” described. So in the face functions shown above, I,J,K would normally be the element local corner number associated with local corner numbers 1,2,3 of the face. Instead we use I =gno(1), J = gno(2), and K = gno(3). In the example above, for the face on the right hand element shared by the two elements in the exploded view, gno(1) = 3, gno(2) = 1, and gno(3) = 2. Note that the discussion here is “1-based” but the code is in c++ so in the code is is “0-based”, e.g. gno(1) is really gno in the code.
This may all seem messy but it is only necessary if you dive deep into the bowels of the basis function routines, such as TetBasisFuncs is file basis.cpp in SRLibSimple. To actually use the library, two setup actions need to be performed:
Create Global Edges For all elements in your code, create corresponding SRelements. This is done by first allocating space for them:
where nel is the number of elements. Then loop over your elements, and, for each, call
model.CreateElem(id, userid, nnodes, nodes, mat)
where nnodes is the number of corner and midnodes in the element, e,g, 10 for a quadratic tet, nodes is the vector of node numbers for the element (corners first followed by midedges), and “mat” contains the element material properties. This are all be described in a programmers guide I uploaded here. The global edges automatically get created on the fly, with the “direction” flag assigned.
Create Global Faces:
After nodes and elements have been defined, just call
This creates global faces for all element faces in the model, and assigns the global node order variable. It also determines which elements own each global face, and which faces are boundary faces (only have one owner).
Subsequently the basis function routines can be called as needs and continuity will be assured.
Homework problem: The volume coordinates for a tetrahedron have the property that they are 1 at a single node, and 0 at the opposite face. For example L1 is 1 at node 1, and 0 on face 2,3,4. Given that property, explain why the blending works properly for the edge, face, and volume functions: For edge functions, the functions are nonzero on the edge, and zero on all other edges and all faces that don’t touch the edge. For face functions, the functions are nonzero on the face, 0 on all other faces. For volume functions, the functions are only nonzero in the interior of the volume.
- Szabo, B, and Babuska, I, Finite Element Analysis, Wiley, NY, 1991