Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Biomech. Author manuscript; available in PMC 2017 March 21.
Published in final edited form as:
PMCID: PMC4801678

Finite Element Simulation of Articular Contact Mechanics with Quadratic Tetrahedral Elements


Although it is easier to generate finite element discretizations with tetrahedral elements, trilinear hexahedral (HEX8) elements are more often used in simulations of articular contact mechanics. This is due to numerical shortcomings of linear tetrahedral (TET4) elements, limited availability of quadratic tetrahedron elements in combination with effective contact algorithms, and the perceived increased computational expense of quadratic finite elements. In this study we implemented both ten-node (TET10) and fifteen-node (TET15) quadratic tetrahedral elements in FEBio ( and compared their accuracy, robustness in terms of convergence behavior and computational cost for simulations relevant to articular contact mechanics. Suitable volume integration and surface integration rules were determined by comparing the results of several benchmark contact problems. The results demonstrated that the surface integration rule used to evaluate the contact integrals for quadratic elements affected both convergence behavior and accuracy of predicted stresses. The computational expense and robustness of both quadratic tetrahedral formulations compared favorably to the HEX8 models. Of note, the TET15 element demonstrated superior convergence behavior and lower computational cost than both the TET10 and HEX8 elements for meshes with similar numbers of degrees of freedom in the contact problems that we examined. Finally, the excellent accuracy and relative efficiency of these quadratic tetrahedral elements was illustrated by comparing their predictions with those for a HEX8 mesh for simulation of articular contact in a fully validated model of the hip. These results demonstrate that TET10 and TET15 elements provide viable alternatives to HEX8 elements for simulation of articular contact mechanics.

Keywords: Finite Element, Biomechanics, Articular Cartilage, Joint, Tetrahedron


Advances in imaging and computational methods make it possible to create and analyze detailed subject-specific models of biomechanical structures from high resolution image data. In this paper we focus on subject-specific finite element (FE) analysis of articular joint contact mechanics. Both generic and subject-specific models of articular contact have been developed and validated to gain insight into load transfer, cartilage mechanics and the etiology of osteoarthritis in the knee (Khoshgoftar et al., 2015; Luczkiewicz et al., 2015), hip (Harris et al., 2012; Henak et al., 2014) ankle (Anderson et al., 2007; Kern and Anderson, 2015) and spine (Dreischarf et al., 2014; Von Forell et al., 2015), among other joints. Despite the progress that has been made in modeling subject-specific joint contact mechanics, many challenges still remain. The articular cartilage of most joints in the human body has complex geometry, undergoes large deformations, is subjected to large compressive loads, and is often thin compared to the surrounding anatomical support. These challenges make it difficult to obtain accurate, validated computational models of articular contact mechanics.

The element type used to discretize the articular geometry is one of the most important choices that affects accuracy and robustness in simulations of articular contact. Linear tetrahedral (TET4) elements are often used due to the ease and robustness of performing automatic meshing, and local and adaptive refinement with tetrahedral elements (Hubsch et al., 1995; Johnson and MacLeod, 1998; Prakash and Ethier, 2001; Spilker et al., 1992) (Delaunay, 1934; Lo, 1991a, b; Lohner, 1996; McErlain et al., 2011; Shephard and Georges, 1991; Wrazidlo et al., 1991). There are several examples in the recent literature that have used TET4 elements to discretize articular cartilage (Das Neves Borges et al., 2014; Johnson et al., 2014; McErlain et al., 2011). However, the TET4 element has several well-known numerical issues. First, TET4 elements can only represent a constant strain state, which necessitates a very fine discretization, often requiring long solution times. Second, TET4 elements lock for nearly incompressible materials as well as under bending deformations (Hughes, 2000), which further reduces their accuracy. Because of these issues, trilinear hexahedral elements (HEX8) have seen much wider use in joint contact analyses, despite the fact that creating hexahedral meshes for complex geometries can be challenging and time consuming. Alternative formulations for TET4 elements have been designed to circumvent their problems (e.g., (Gee et al., 2009; Puso and Solberg, 2006) based on nodal averaging of the deformation gradient). Although they reduce locking, they have other problems such us spurious deformation modes (Maas, 2011), which make them inaccurate for contact analysis.

Quadratic tetrahedral elements are an attractive alternative to TET4 elements. They maintain the advantages of tetrahedral mesh generation and they can represent curved boundaries more accurately than HEX8 elements since their edges and faces can deform. Although quadratic tetrahedral elements have been investigated as alternatives to HEX8 elements (e.g., (Cifuentes and Kalbag, 1992; Tadepalli et al., 2011; Weingarten, 1994)), none of these studies investigated their application for simulation of articular contact analyses. The 10-node tetrahedron (TET10) has seen some limited use for contact mechanics (Bunbar et al., 2001; Hao et al., 2011; Tadepalli et al., 2011; Wan et al., 2013; Yang and Spilker, 2006). To our knowledge the 15-node tetrahedron (TET15) has never been used in computational biomechanics and it has seen very limited use in nonlinear computational solid mechanics at all (Danielson, 2014), although it has been used in the fluid mechanics community (Bertrand et al., 1992).

The objectives of this study were to determine the efficacy in terms of accuracy of the recovered stresses, robustness in terms of the convergence behavior, and computational expense of the TET10 and TET15 elements compared to the HEX8 element in the context of articular contact mechanics. First, the effects of different integration rules on the stress predictions and computational cost were investigated and used to determine suitable integration rules for both the TET10 and TET15 elements. Then using these integration rules, the accuracy and computational cost of both elements were compared to the HEX8 element for several benchmark contact problems. For one of these problems, we examined the results obtained using TET4 elements to contrast their performance with the quadratic elements. Finally, we compared predictions of contact stresses using HEX8, TET10, and TET15 elements for a validated model of hip contact mechanics (Henak et al., 2014).


Element formulation and numerical integration

All element formulations in this research were implemented in the FEBio software suite (, which uses an implicit, Newton-based method to solve the nonlinear FE equations of solid mechanics (Maas et al., 2012a). The TET10 element has 10 nodes: 4 corner nodes and 6 nodes located at the midpoint of the edges (Figure 1). Due to the quadratic shape functions, the facets and edges of this element can distort and therefore the element behavior is “softer” than the TET4 element. The TET15 element adds one more node at the center of each facet, and one in the center of the tetrahedron (Figure 1). Although the TET15 element has more nodes than the TET10, it is still a quadratic element since the highest order of complete polynomial that can be represented by the shape functions is second order. The TET15 element can represent some forms of quadratic strains, whereas the TET10 element can only represent linear strains.

Figure 1
Schematic of the node topology for two quadratic tetrahedral elements that were examined in this study. Left, 10-node quadratic tetrahedron (TET10). Right, 15-node quadratic tetrahedron (TET15). Closed circles represent corner nodes, open circles represent ...

In a FE formulation, the discretized form of the equilibrium equations requires the use of appropriate numerical integration schemes. For second-order elements, the integration rule should have at least second-order accuracy, at least for linear analyses. The ideal integration rule is less clear for large deformation nonlinear analyses. For this reason, we implemented and compared several volumetric integration rules. We denote volume integration rules as V(n) where (n) indicates the number of integration points. Similarly, the notation S(n) is used to denote surface integration rules. For the TET10 element, both 4-point (V4) and 8-point (V8) Gauss integration rules were implemented (Abramowitz and Stegun, 1964). For the TET15 element, 11-point (V11) and 15-point (V15) Gauss integration rules were implemented (Keast, 1986). All of these rules are at least second-order accurate (for linear analyses) and are symmetric, i.e. the integration points are distributed in a symmetrical spatial pattern.

For contact enforcement, a surface integration rule is required to integrate the traction forces over the discrete surface, represented by facets of the finite elements. The facets of the TET10 element are 6-node quadratic triangles. Two surface integration rules were implemented and compared: 3-point Gauss (S3), and 7-point Gauss (S7) (Figure 2). For the TET15 element, a 7-node quadratic facet, and the same S3 and S7 integration rules were implemented. Again, these rules have at least second-order accuracy and are symmetric (Abramowitz and Stegun, 1964).

Figure 2
Schematic of the surface integration rules that were implemented and tested for evaluation of the contact integrals for the quadratic triangle facets composing the surface of the tetrahedral elements. Left panel: S3 integration rule. Right panel: S7 integration ...

The mean dilatation formulation was used for the HEX8 models (Simo and Taylor, 1991). This formulation is known to perform well for nearly-incompressible materials for which the standard displacement-based formulation, using full integration, has a tendency to lock.

All analyses were performed using the quasi-Newton solver in FEBio (Maas et al., 2012b), which is based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method (Matthies and Strang, 1979). This nonlinear solver begins with a full formation and factorization of the stiffness matrix and then proceeds with a user-defined number of BFGS updates, which involve computation of the right-hand-side (RHS) vector. For two of the problems below, we report both the number of stiffness matrix reformations and RHS evaluations as metrics of nonlinear convergence and computational effort. All contact analyses used the “sliding-tension-compression” contact algorithm in FEBio (Maas et al., 2012b). This algorithm implements a facet-on-facet, frictionless sliding contact where the contacting surfaces can separate but not penetrate (Laursen, 2002). The augmented Lagrangian method was used to enforce the contact constraint to a user-defined tolerance (Laursen and Maker, 1995).

Stress recovery

During the FE solution process, stresses are typically evaluated at the integration points. Since stresses are calculated using the shape function derivatives, they are usually discontinuous across element boundaries. For visualization and for further postprocessing analysis (e.g. a-posteriori stress error estimates (Zienkiewicz and Zhu, 1992)), a more accurate, continuous stress field is required. This is accomplished using a stress-recovery algorithm that projects the stresses from the integration points to the nodes. For TET4 elements, a weighed nodal average of the surrounding element values is usually sufficient, but for quadratic elements this does not produce good results. We implemented and used the Superconvergent Patch Recovery (SPR) method to recover stresses for the quadratic elements (Zienkiewicz and Zhu, 1992). This method incorporates a least-squares approach to fit a polynomial to the integration point values of an element patch surrounding each node. The polynomial is then evaluated at the nodes to determine the nodal value. Another important reason for investigating different integration rules is to determine the minimal integration rule that recovers the stresses accurately. Since the SPR method uses least-squares, each node must be surrounded by a minimal number of integration points to solve the least-squares system. This could pose problems at the boundary of the mesh and in thin structures. For the TET10 element, we used a linear least-squares fit since there are often not enough integration points (especially for the V4 integration rule) to fit a complete quadratic polynomial. For the TET15 element, a complete quadratic polynomial was used.

Plane strain contact

Different combinations of surface and volume integration rules were investigated using a plane strain contact model (Figure 3). This model is representative of articular contact mechanics in that it models the contact between curved, thin, incongruent, deformable materials. The model consisted of two parts. The top part had dimensions of 1.5 mm (height)×20 mm (length)×1 mm (thickness). The bottom part had the same height and thickness but was curved downward with a radius of curvature of 58 mm. The bottom of the curved layer was fixed while the top of the flat layer was displaced downward. Since FEBio is designed for 3D analysis, plane strain boundary conditions were created by preventing the nodes from displacing in the out-of-plane direction. The materials of both layers were represented by an uncoupled Mooney-Rivlin constitutive model with strain energy:

Figure 3
Schematic of the plane strain contact problem. The bottom surface of the curved layer is held fixed, while the top layer undergoes a prescribed downward displacement. Dashed area shows the region illustrated in the contour plots in Figure 5.

Here, Ĩ1, Ĩ2 are the first and second invariants of the deviatoric right Cauchy-Green deformation tensor, and J is the determinant of the deformation gradient. This material formulation and the material coefficients below were chosen to represent articular and femoral articular cartilage and coincide with the properties of a validated hip joint model (discussed below). The material parameters were c1 = 2.98 MPa, c2 = 0 MPa, k = 594 MPa, for the flat layer; and c1 = 1.86 MPa, c2 = 0 MPa, k = 371 MPa for the curved layer. To compare the integration rules, we examined isovalue fringe plots of the 3rd principal stress and searched for erroneous patterns in the isovalues.

Compression and sliding of a rigid cylinder on a deformable block

To assess the performance of the quadratic tetrahedral elements in contact problems with large amounts of sliding, a benchmark problem from the FEBio Verification Suite was used (co07.feb). This is a revelant test problem since many joints (e.g. knee, shoulder, hip) undergo large amounts of sliding during articular contact. A rigid cylinder (outer radius = 1 mm) first indents a rectangular block (length = 5 mm, width = 1.9 mm, height = 1 mm) and then slides across the surface of the block (Figure 4). The cylinder was modeled as a rigid body using quadratic 20-node hexahedral elements. Mid-side nodes of the elements were positioned so that the surface of the cylinder was smooth to avoid jumps in the contact stiffness due to discontinuities in the surface normal. Again, as before, the material formulation was chosen to represent articular cartilage. The rectangular block was represented using the Mooney-Rivlin model above (c1 = 2.98 MPa, c2 = 0 MPa, k = 594 MPa) and was discretized with TET4, TET10 and TET15 elements. A model using HEX8 elements was created and used as a reference. A mesh convergence study was conducted to allow fair assessment of the computational costs of the four element types (Table 1). Convergence was measured in terms of average reaction force on the rigid indenter over the entire horizontal displacement of the indenter. Fringe plots of the 3rd principal stress at the midpoint of the analysis were used to compare the predictions of the models.

Figure 4
Geometry and loading path of the rigid cylinder on the deformable box. First, the cylinder was displaced vertically. After reaching maximal vertical displacement, a horizontal displacement was applied while the vertical position was maintained. From top ...
Table 1
Results of the mesh convergence study for the contact sliding problem showing the number of nodes and elements, the number of time steps to complete the analysis, the number of matrix reformations (Refs) and right-hand side (RHS) evaluations, and runtime ...

Articular contact in the human hip

To test the effectiveness of quadratic tetrahedral elements for production size articular contact problems involving complex geometry, the articular layers in a model of the hip from a previous validation study (Henak et al., 2014) were discretized with tetrahedral elements (Figure 9). In this previous study, articular layers were discretized with HEX8 elements and after a mesh convergence study it was determined that six elements through the thickness produced converged results. The thickness varied across the articular surfaces and ranged from 0.2 mm to 3.7 mm on the femoral side and from 0.5 mm to 4.5 mm on the acetabular side. The HEX8 model was validated with experimental measurements of contact stress. Cartilage layers were modeled using a Mooney-Rivlin material (acetabular cartilage c1=2.98 MPa, c2=0 MPa, k=594 MPa, femoral cartilage c1=1.86 MPa, c2=0 MPa, k=371 MPa). For this study, two additional models were created using TET10 and TET15 elements. The mesh resolution was chosen such that the number of nodes (and thus the degrees of freedom) roughly matched the HEX8 model (Table 2). For the tetrahedral models the number of elements through the thickness varied from one to four depending on the cartilage thickness. In the areas where contact was occurring there were about two to three elements through the thickness on average. Loading conditions were applied that simulated walking heel strike (Bergmann et al., 2001). A prescribed load of 1650 N was applied to the femur while keeping the pelvis fixed. The 3rd principal stress, an indicator of contact stress, and contact area were compared.

Figure 9
A) Finite element model of the human hip, including the acetabulum, proximal femur and articular layers on the acetabulum and femoral head. Red square shows approximate area corresponding to contour plots. Blue line shows approximate location of inset, ...
Table 2
Result of the hip model study where the results of a validated HEX8 model were compared with the predictions of a TET10 and TET15 model. The table lists the number of nodes and elements, peak contact stress, the contact area, the number of right-hand-side ...


Plane strain contact

The surface integration rule had a large effect on predictions of 3rd principal stress (Figure 5). The S3 surface integration rule predicted overall larger variations in stress near the contact interface than the S7 rule (compare rows 2 and 3 in Figure 5). We concluded that the S7 surface integration rule is the best choice of the two rules. In contrast, there was little effect of volume integration rule on predicted stresses (compare rows 3 and 4 in Figure 5). Indeed, comparing for instance the combination S3/V4 with S3/V8 for the TET10 element (Figure 5, left column), the results are nearly indistinguishable. A noticeable exception is the combination S3/V11 and S3/V15 for the TET15 element where the stresses do show a significant difference. We found that the combination S3/V11 ran very poorly and we suspect that this combination is unstable and should be avoided. Although the other results show little variation for the different volume integration rules, we did observe improved convergence behavior for the higher-order volume integration rules. For this reason the remaining simulations were done with the S7/V8 rule for the TET10 element and the S7/V15 for the TET15 element.

Figure 5
Results for the plane strain contact problem. The isovalue fringe plots show 3rd principal stress for various combinations of surface and volume integration rules. Left column: results for TET10 element. Right column: results for the TET15 element. The ...

Compression and sliding of a rigid cylinder on a deformable block

Overall reaction force on the cylinder provided an integrated measure of mesh convergence (Figure 6). The quadratic tetrahedral element formulations demonstrated excellent agreement above 10,000 nodes with the HEX8 results. With the exception of the coarsest meshes, all values are within one percent of each other (Table 1, far right column). This is in stark contrast with the TET4 results, which did not appear converged even at the finest resolution and predicted much larger reaction forces.

Figure 6
Average reaction force on the rigid indenter while the indenter slides horizontally across the block, maintaining maximal vertical displacement. Average reaction forces are shown as a function of the number of nodes for a TET4 mesh, a TET10 mesh, a TET15 ...

Inspecting the reaction force as a function of horizontal displacement can give an indication of the sensitivity of the reaction force to the mesh distribution. At the lowest mesh resolution all formulations showed fluctuations in the reaction forces, however at the finest resolution the reaction force varied smoothly (Figure 7). The difference between the minimum and maximum value was about 3% for each element formulation. The differences between the HEX8, TET10, and TET15 element formulations at a given displacement were even smaller and were less than 1%. Note again that the TET4 model predicted signifcantly larger reaction forces than the other three formulations.

Figure 7
Variation of the reaction force on the rigid indenter as a function of horizontal displacement for the finest mesh resolution of HEX8, TET4, TET10, and TET15 models. The variation is relatively small and smooth throughout the entire displacement and is ...

Converence behavior was comparable between the four element formulations. The number of stiffness matrix reformations and right-hand side evaluations were similar. The run times for the HEX8 were similar, though consistently higher, than those for the TET4 and TET10 formulations. The TET15 models ran sigfnicantly faster. The predicted 3rd principal stress compared well between the quadratic tetrahedral and hexahedral models at the finest mesh resolutions (Figure 8, right column). For the coarsest mesh resolutions, the TET15 and TET10 models were similar to the HEX8 results, though variations due to the irregularity of the mesh were noticeable (Figure 8). The TET4 results were significantly higher and showed mesh dependent variations at all mesh resolutions.

Figure 8
3rd principal stress fringe plots for two different mesh densities at the half-way point of horizontal displacement. From left to right: coarsest mesh, finest mesh density. From top to bottom: HEX8, TET4, TET10, and TET15 mesh. Stress predictions were ...

Articular Contact in the Hip

Based on the results of the previous two models, the S7 surface integration rule was used to evaluate the contact integrals for both tetrahedral elements in this model. The V8 and V15 volume integration rules were used for the TET10 and TET15 elements, respectively. There was good qualitative agreement of the contact stress (measured by the 3rd principal stress) between the three models (Figure 9). Quantitative comparisons confirmed this observation. For both the TET10 and TET15 models, the peak contact stress and the contact area on the acetabular cartilage compared well to the HEX8 results (Table 2).

In terms of performance, measured by the convergence statistics as well as overall runtime, all three models ran comparably, although the TET10 model did run slightly longer than both the HEX8 and TET15 models (Table 2).


The results of this study demonstrate that quadratic tetrahedral elements are comparable both in computational cost and accuracy to the “gold standard” HEX8 elements for articular contact analysis. In some cases the TET15 models ran faster than the other models with comparable numbers of degrees of freedom, while in other cases they ran only slightly slower. It may seem surprising that these elements with more degrees of freedom are competitive in terms of computational cost, and in fact are faster than the TET4 model (e.g. the sliding contact example). Although the evaluation of the element quantities (i.e. element residual vector and stiffness matrix) is faster for low-order elements, the overall computational time may not follow the same pattern, as the examples clearly demonstrated. This can be explained by recognizing that many operations in the FE method (evaluation of global quantities, matrix assembly, etc.) are more proportional to the number of integration points than the number of degrees of freedom. This is especially true for contact problems, where the projection algorithm that determines the location of contact between surfaces can constitute a significant amount of the solution time. For a given number of nodes, a TET15 model will always have fewer elements than a similar TET10 or HEX8 and, although it uses more integration points than the TET10 or HEX8, the net sum often appears to be a reduction in the number of integration points. This may in part explain why the TET15 models ran faster in many cases.

The results of the plane strain contact analyses were very sensitive to the surface integration rule, and the S7 integration rule produced good results for the quadratic tetrahedrons. A possible explanation for this can be found by realizing that the accuracy of (Gaussian) integration rules is usually related to the ability to integrate polynomials exactly. For contact problems with non-conforming and possibly distorted surfaces, polynomials may not fit the surface tractions well and thus higher order integration rules are necessary to integrate the contact tractions accurately.

In contrast, the volume integration rule had little effect on the accuracy of the recovered stresses. All volume integration rules that were investigated are sufficient to obtain accurate results in linear analysis and it appears that, for the problems considered, the integration rules were also sufficiently accurate for nonlinear contact analysis. We did observe that the S3 surface integration rule in combination with the V11 volume integration rule for TET15 elements produced poor convergence as well as inaccurate results. We suspect that this combination does not lead to a stable element formulation. Despite the relatively small effect of the volume integration rule on the stress predictions, it may be prudent to consider the higher-order rules. The SPR recovery method, as with any least-squares method, requires a sufficient number of sample points to make the problem well-posed. In the present context this amounts to having a sufficient number of integration points surrounding each node. This may be challenging near boundaries, especially for articular contact problems, since the cartilage layers generally tend to be relatively thin. The result is that the V4 volume integration rule for TET10 prevents the use of a truly quadratic fit in the SPR recovery method. Although in this study we used a linear fit for the TET10 element and thus the V4 rule could have sufficed, there may be situations where the V8 rule is necessary to achieve accurate results. For the same reasons, for the TET15 elements, which contain additional monomials in their shape functions and thus require even more integration points that the TET10, it is prudent to consider the higher-order volume integration rule.

In this study we used “traditional” facet-on-facet based contact formulations (Laursen, 2002). In recent years, it has been demonstrated that mortar-based contact formulations have several advantages over traditional contact formulations (Puso et al., 2008). In essence, they minimize the error introduced by the non-conforming interface in either the displacements or the contact tractions. Yet, traditional formulations are widely used in many open-source and commercial finite element packages and thus evaluation of quadratic tetrahedral elements with traditional contact formulations has significant practical merit. These traditional formulations are known to give satisfactory results and, when used with an augmented Lagrangian method offer a good compromise between computational cost and accuracy. In addition, we anticipate that even for mortar-contact formulations the use of quadratic tetrahedral formulations will offer a distinctive advantage, foremost due their ability to represent arbitrarily curved surfaces more accurately than linear triangular elements. Thus, we expect that the findings in this study remain valid, regardless of the contact formulation, though further investigations are necessary to confirm this.

In this study we only looked at the effects of frictionless contact. This is not seen as a strong limitation since friction is usually not considered in articular contact. That said, friction is a feature that is still under development in FEBio’s contact formulations. Similarly, the effect of fluid motion in the tissue is of considerable interest in articular contact and was not considered here. The incorporation of fluid motion requires a biphasic analysis and the performance of quadratic tetrahedral formulations in this context is an important area for future investigation.

In the sliding contact case study we compared the quadratic tetrahedron results with linear tetrahedron results and demonstrated that, given the same number of nodes, the quadratic formulations, and especially the TET15 element, performed much better in all areas (i.e. accuracy, convergence, run-times) than the linear tetrahedron. In a way, this can be seen as yet another confirmation of the well-known mantra in finite element analysis that so-called p-refinement (i.e. increasing the polynomial order) is often a better refinement strategy than h-refinement (i.e. reducing the element size). Simply moving from a linear to a quadratic formulation resulted in vast improvements. In addition, the decreased run-times of higher-order tetrahedral formulations (for a given number of nodes), although perhaps surprising, is seen as an additional advantage of these quadratic tetrahedral formulations.

In conclusion, quadratic tetrahedral formulations offer an excellent alternative to HEX8 elements for simulation of articular contact. The TET15 element in particular appears to be very well suited for articular contact analysis.


Financial support from NIH grants #R01AR053344 and #R01GM083925 is gratefully acknowledged.


Conflict of Interest Statement

The authors certify that they have no conflict of interest in relation to the submitted manuscript.

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


  • Abramowitz M, Stegun IA. Handbook of mathematical functions. 1964.
  • Anderson DD, Goldsworthy JK, Li W, James Rudert M, Tochigi Y, Brown TD. Physical validation of a patient-specific contact finite element model of the ankle. J Biomech. 2007;40:1662–1669. [PMC free article] [PubMed]
  • Bergmann G, Deuretzbacher G, Heller M, Graichen F, Rohlmann A, Strauss J, Duda GN. Hip contact forces and gait patterns from routine activities. Journal Of Biomechanics. 2001;34:859–871. [PubMed]
  • Bertrand FH, Gadbois MR, Tanguy PA. Tetrahedral elements for fluid flow. International Journal for Numerical Methods in Engineering. 1992;33:1251–1267.
  • Bunbar WL, Un K, Donzelli PS, Spilker RL. An evaluation of three-dimensional diarthrodial joint contact using penetration data and the finite element method. Journal of Biomechanical Engineering. 2001;123:333–340. [PubMed]
  • Cifuentes AO, Kalbag A. A performance study of tetrahedral and hexahedral elements in 3-D finite element structural analysis. Finite Elements in Analysis and Design. 1992;12:313–318.
  • Danielson KT. Fifteen node tetrahedral elements for explicit methods in nonlinear solid dynamics. Computer methods in applied mechanics and engineering. 2014;272:160–180.
  • Das Neves Borges P, Forte AE, Vincent TL, Dini D, Marenzana M. Rapid, automated imaging of mouse articular cartilage by microCT for early detection of osteoarthritis and finite element modelling of joint mechanics. Osteoarthritis Cartilage. 2014;22:1419–1428. [PMC free article] [PubMed]
  • Delaunay BNSlSVIANS. VII Seria Otdelenie Matematicheskii i Estestvennyka Nauk. 1934;7:793–800. “Sur la Sphere” Vide” Izvestia Akademia Nauk SSSR, VII Seria, Otdelenie Matematicheskii i Estestvennyka Nauk 7 793 800.
  • Dreischarf M, Zander T, Shirazi-Adl A, Puttliz CM, Adam CJ, Chen CS, Goel VK, Kiapour A, Kim YH, Labus KM, Little JP, Park WM, Wang YH, Wilke HJ, Rohlmann A, Schimdt H. Comparison of eight published static finite element models of the intact lumbar spine: Predictive power of models improves when combined together. Journal of Biomechanics. 2014;47:1757–1766. [PubMed]
  • Gee MW, Dohrmann CR, Key SW, Wall WA. A uniform nodal strain tetrahedron with isochoric stabilization. International Journal for Numerical Methods in Engineering. 2009;78:429–443.
  • Hao Z, Wan C, Gao X, Ji T. The effect of boundary condition on the biomechanics of a human pelvic joint under an axial compressive load: A three-dimensional finite element model. Journal of Biomechanical Engineering. 2011:133. [PubMed]
  • Harris MD, Anderson AE, Henak CR, Ellis BJ, Peters CL, Weiss JA. Finite element prediction of cartilage contact stresses in normal human hips. J Orthop Res. 2012;30:1133–1139. [PMC free article] [PubMed]
  • Henak CL, Kapron AL, Anderson AE, Ellis BJ, Weiss JA. Specimen-specific predictions of contact stress under physiological loading in the human hip: validation and sensitivity studies. Biomech Model Mechanobiol. 2014;13:387–400. [PMC free article] [PubMed]
  • Hubsch PF, Middleton J, Meroi EA, Natali AN. Adaptive finite-element approach for analysis of bone/prosthesis interaction. Med Biol Eng Comput. 1995;33:33–37. [PubMed]
  • Hughes TJR. The Finite Element Method. Dover; 2000.
  • Johnson CR, MacLeod RS. Adaptive local regularization methods for the inverse ECG problem. Prog Biophys Mol Biol. 1998;69:405–423. [PubMed]
  • Johnson JE, Lee P, McIff TE, Bruce Toby E, Fischer KJ. Computationally Efficient Magnetic Resonance Imaging Based Surface Contact Modeling as a Tool to Evaluate Joint Injuries and Outcomes of Surgical Interventions Compared to Finite Element Modeling. Journal of Biomechanical Engineering. 2014;136:041002–041002. [PubMed]
  • Keast P. Moderate-degree Tetrahedral Quadrature Formulas. Computer methods in applied mechanics and engineering. 1986;55:339–348.
  • Kern AM, Anderson DD. Expedited patient-specific assessment of contact stress exposure in the ankle joint following definitive articular fracture reduction. J Biomech 2015 [PMC free article] [PubMed]
  • Khoshgoftar M, Vrancken AC, van Tienen TG, Buma P, Janssen D, Verdonschot N. The sensitivity of cartilage contact pressures in the knee joint to the size and shape of an anatomically shaped meniscal implant. J Biomech. 2015;48:1427–1435. [PubMed]
  • Laursen TA. Computational Contact and Impact Mechanics. Springer-Verlag; 2002.
  • Laursen TA, Maker BN. Augmented Lagrangian Quasi-Newton Solver for Constrained Nonlinear Finite Element Applications. International Journal for Numerical Methods in Engineering. 1995;38:3571–3590.
  • Lo SH. Volume Discretization into Tetrahedra-I. Verification and Orientation of Boundary Surfaces. Computers and Structures. 1991a;39:493–500.
  • Lo SH. Volume Discretization into Tetrahedra - II. 3D Triangulation by Advancing Front Approach. Computers and Structures. 1991b;39:501–511.
  • Lohner R. Progress in Grid Generation via the Advancing Front Technique. Engineering with Computers. 1996;12:186–210.
  • Luczkiewicz P, Daszkiewicz K, Witkowski W, Chroscielewski J, Zarzycki W. Influence of meniscus shape in the cross sectional plane on the knee contact mechanics. J Biomech. 2015;48:1356–1363. [PubMed]
  • Maas SA. Implementation and verification of a nodally-integrated tetrahedral element in FEBio. SCI Technical Report 2011
  • Maas SA, Ellis BJ, Ateshian GA, Weiss JA. FEBio: Finite Elements for Biomechanics. JBME. 2012a:134. [PubMed]
  • Maas SA, Rawlins DS, Weiss JA, Ateshian GA. FEBio User’s Manual 1.5 2012b
  • Matthies H, Strang G. The solution of nonlinear finite element equations. Int J Numer Methods Eng. 1979;14:1613–1626.
  • McErlain DD, Milner JS, Ivanov TG, Jencikova-Celerin L, Pollmann SI, Holdsworth DW. Subchondral cysts create increased intra-osseous stress in early knee OA: A finite element analysis using simulated lesions. Bone. 2011;48:639–646. [PubMed]
  • Prakash S, Ethier CR. Requirements for mesh resolution in 3D computational hemodynamics. J Biomech Eng. 2001;123:134–144. [PubMed]
  • Puso M, Solberg J. A stabilized nodally integrated tetrahedral. International Journal for Numerical Methods in Engineering. 2006;67:841–867.
  • Puso MA, Laursen TA, Solberg J. A segment-to-segment mortar contact method for quadratic elements and large deformations. Computer methods in applied mechanics and engineering. 2008;197:555–566.
  • Shephard MS, Georges MK. Three-Dimensional Mesh Generation by Finite Octree Technique. International Journal for Numerical Methods in Engineering. 1991;32:709–749.
  • Simo JC, Taylor RL. Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms. Computer methods in applied mechanics and engineering. 1991;85:273–310.
  • Spilker RL, de Almeida ES, Donzelli PS. Finite element methods for the biomechanics of soft hydrated tissues: nonlinear analysis and adaptive control of meshes. Crit Rev Biomed Eng. 1992;20:279–313. [PubMed]
  • Tadepalli SC, Erdemir A, Cavanagh PR. Comparison of hexahedral and tetrahedral elements in finite element analysis of the foot and footwear. Journal of Biomechanics. 2011;44:2337–2343. [PubMed]
  • Von Forell G, Stephens T, Samartzis D, Bowden A. Low Back Pain: A Biomechanical Rationale Based on “Patterns” of Disc Degeneration. Spine. 2015;40:1165–1172. [PubMed]
  • Wan C, Hoa Z, Wen S. The effect of the variation in ACL constitutive model on joint kinematics and biomechanics under different loads: a finite element study. Journal of Biomechanical Engineering. 2013;135:041002. [PubMed]
  • Weingarten VI. The controversy over hex or tet meshing. Machine Design. 1994;66:74–78.
  • Wrazidlo W, Brambs HJ, Lederer W, Schneider S, Geiger B, Fischer C. An alternative method of three-dimensional reconstruction from two- dimensional CT and MR data sets. Eur J Radiol. 1991;12:11–16. [PubMed]
  • Yang T, Spilker RL. A Lagrange multiplier mixed finite element formulation for three-dimensional contact of biphasic tissues. Journal of Biomechanical Engineering. 2006;129:457–471. [PubMed]
  • Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique. International Journal for Numerical Methods in Engineering. 1992;33:1331–1634.