The stiffness matrix for a conventional isoparametric finite element takes the form:
Kel = ∫BTCBdV
where C is the element stiffness or constitutive matrix, and B expresses the strain-displacement relations. If the strain is stored as a 1D vector e = [ex,e
y,ez,γxy,yxz,γyz]. and displacements as u = [ux,uy,uz], then
e = [B]u.
The strain-displacement matrix B has 3×6 submatrices for each basis function hj:
Evaluating B is the only difference between a conventional element routine and an element with p-adaptive functions. In a conventional routine, the shape functions are used for the mapping, and the same functions are used for the basis functions. In a p-adaptive element, different functions are used for mapping and displacement. In stressRefine, quadratic serendipity functions are used for the mapping, and higher order polynomials for the displacements.
where r,s,t are the natural coordinates in an element and J-1 is the inverse of the mapping from r,s,t at each integration point. The integral over the volume is also converted into an integral in natural coordinates using the determinant of the mapping |J|. J is computed from the shape functions. It is invertible to J-1 unless the element is invalid (e.g. too highly distorted).
Examining the code.
A simple implementation of a p-adaptive stiffness matrix is in the routine CalculateStiffnessMatrix in SRelement.cpp in the stressRefine library SRlibSimple.
int nfun = globalFunctionNumbers.GetNum();
looks up the number of displacement basis functions in the element, which is calculated from the polynomial order of each edge of the element.
The integration points for the element are determined with the call to model.math.FillGaussPoints. The stressRefine library uses degenerated brick Gauss quadrature for tets and wedges. This could be made more efficient by using the triangular points developed by Cowper  for the r,s quadrature in wedges and the tetrahedral points developed by JinYun . These do not go up to the higher polynomial orders needed by stressRefine, but could be used when the polynomial order is low enough, switching to the degenerated Gauss quadrature when needed.
There is a quadrature loop in the element matrices
for (gp = 0; gp < nint; gp++)
Inside that loop, the natural coordinates and the quadrature weight are determined in model.math.GetGP3d, and J-1 and |J| are calculated in FillMapping. An error is raised if |J| is too small. However, before the element routines are calculated, stressRefine tests the mapping of each element and attempts to recover by partially flattening curved elements as discussed previously.
The derivatives of the basis functions ∂hj/∂r, etc are calculated at each integration point with the call FillBasisFuncs.
After this the x,y and z derivatives of the basis functions can be calculated.
BTC is calculated in fillBTC, which accounts for the zeros in B and C.
The multiplication (BTC) times B is then calculated in FillKel33GenAnisoRowWise.
This returns a 3×3 submatrix kel33 of the element matrix which is stored in the appropriate location in the symmetcially-stored element stiffness matrix.
Converting an existing conventional element stiffness routine
Any stiffness routine will have a loop over the integration points, and loops for the rows and columns over the number of basis functions in the element.
A few modifications need to be made. This discussion assumes the setup functions for stressRefine have been called, and for each conventional element for which p-adaptivity will be used, a corresponding SRelement has been created, This was discussed in lesson 2. Then CountElementFunctions can be called in the library function in SRbasis to determine the number of element functions. FillGaussPoints in model.math will calculate and store the integration points and return the number. The numerical integration in the element must be modified to use the increased number appropriate to the number of basis functions corresponding to the polynomial order of the element.
The derivatives of the higher-order polynomial basis functions ∂hj/∂r, etc must be used instead of the shape functions. These can be calculated directly using ElementBasisFuncs in SRbasis. Everything else in the element routine is unchanged.
The stresses in an element can be directly computed from
σ = [D]ε = = [D][B]u once the displacements u are known.
This is computed the same way as for a conventional element, except, as for the case with the element stiffness matrix, the displacement gradients in [B] are calculated using the higher-order polynomial basis functions, not the conventional shape functions.
Adapting the polynomial order
Between this and the previous lesson, we can calculate the stiffness matrix for an element given the polynomial orders of each of its edges. In the next lesson we’ll cover adaptivity: after determining the solution for the current polynomial order, errors are calculated, and used to estimate the required polynomial order to achieve the desired accuracy.
Homework problem: Is the element routine CalculateStiffnessMatrix “thread-safe” if we want to calculate the elements in parallel?
- Cowper, G, “Gaussian Quadrature Formulae For Triangles”, Int J Num Meth Engr, 1973
- Jinyun, Yu, “Symmetric Gaussian Quadrature Formulae For Tetrahedral Regions”, Computer Meth. Appl. Mech. Engr, 1984.