Chapter 8

 

Theory and Implementation of the Finite Element Method

 

 

 

8.6 Advanced element formulations $–$ Incompatible modes; reduced integration; and hybrid elements

 

Techniques for interpolating the displacement field within 2D and 3D finite elements were discussed in Section 8.1.9 and 8.1.10.  In addition, methods for evaluating the volume or area integrals in the principle of virtual work were discussed in Section 8.1.11.   These procedures work well for most applications, but there are situations where the simple element formulations can give very inaccurate results.  In this section

1. We illustrate some of the unexpected difficulties that can arise in apparently perfectly well designed finite element solutions to boundary value problems;
2. We describe a few more sophisticated elements that have been developed to solve these problems.

 

We focus in particular on `Locking’ phenomena.  Finite elements are said to `lock’ if they exhibit an unphysically stiff response to deformation.  Locking can occur for many different reasons. The most common causes are (i) the governing equations you are trying to solve are poorly conditioned, which leads to an ill conditioned system of finite element equations; (ii) the element interpolation functions are unable to approximate accurately the strain distribution in the solid, so the solution converges very slowly as the mesh size is reduced; (iii) in certain element formulations (especially beam, plate and shell elements) displacements and their derivatives are interpolated separately. Locking can occur in these elements if the interpolation functions for displacements and their derivatives are not consistent.

 

 

8.6.1 Shear locking and incompatible mode elements

 

Shear locking can be illustrated by attempting to find a finite element solution to the simple boundary value problem illustrated in the picture. Consider a cantilever beam, with length L, height 2a and out-of-plane thickness b, as shown in the figure.  The top and bottom of the beam $x_2=\pm a$ are traction free, the left hand end is subjected to a resultant force P, and the right hand end is clamped.  Assume that b<<a, so that a state of plane stress is developed in the beam.  The analytical solution to this problem is given in Section 5.2.4.

 

The figures below compare this result to a finite element solution, obtained with standard 4 noded linear quadrilateral plane stress elements.  Results are shown for two different ratios of $a/L$, with $P{L}^3/E{a}^{3}b=2.22$ for both cases. 

   

 

For the thick beam, the finite element and exact solutions agree nearly perfectly.   For the thin beam, however, the finite element solution is very poor $–$ even though the mesh resolution is unchanged. The error in the finite element solution occurs because the standard 4 noded quadrilateral elements cannot accurately approximate the strain distribution associated with bending.  The phenomenon is known as `shear locking’ because the element interpolation functions give rise to large, unphysical shear strains in bent elements.  The solution would eventually converge if a large number of elements were added along the length of the beam. The elements would have to be roughly square, which would require about 133 elements along the length of the beam, giving a total mesh size of about 500 elements. 

 

Shear locking is therefore relatively benign, since it can be detected by refining the mesh, and can be avoided by using a sufficiently fine mesh.  However, finite element analysts sometimes cannot resist the temptation to reduce computational cost by using elongated elements, which can introduce errors.

 

Shear locking can also be avoided by using more sophisticated element interpolation functions that can accurately approximate bending. `Incompatible Mode’ elements do this by adding an additional strain distribution to the element.  The elements are called `incompatible’ because the strain is not required to be compatible with the displacement interpolation functions. The approach is conceptually straightforward:

1.      The displacement fields in the element are interpolated using the standard scheme, by setting

$u_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)u_i^a}\delta v_i(x)={\displaystyle \sum _{a=1}^n{N}^a(x)\delta v_i^a}{x}_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)x_i^a}$ 

where $N^a({\xi }_j)$ are the shape functions listed in Sections 8.1.9 or 8.1.10, ${\xi }_i$ are a set of local coordinates in the element, $u_i^a,x_i^a$ denote the displacement values and coordinates of the nodes on the element, and $N_e$ is the number of nodes on the element.

2.      The Jacobian matrix for the interpolation functions, its determinant, and its inverse are defined in the usual way

$\frac{\partial x_i}{\partial {\xi }_j}={\displaystyle \sum _{a=1}^{N_e}\frac{\partial N^a(\xi )}{\partial {\xi }_j}{x}_i^a}$    $J=\det \left(\frac{\partial x_i}{\partial {\xi }_j}\right)$        $\frac{\partial {\xi }_j}{\partial x_i}={\left(\frac{\partial x}{\partial \xi }\right)}_{ji}^{-1}$

3.      The usual expression for displacement gradient in the element is replaced by

$\frac{\partial u_i}{\partial x_j}=\left({\displaystyle \sum _{a=1}^{N_e}\frac{\partial N^a(\xi )}{\partial {\xi }_m}}{u}_i^a+{\displaystyle \sum _{k=1}^p\frac{J(0)}{J(\xi )}}{\alpha }_i^{(k)}{\delta }_{km}{\xi }_k\right)\frac{\partial {\xi }_m}{\partial x_j}$

where p=2 for a 2D problem and p=3 for a 3D problem, ${\alpha }_i^{(k)}$ are a set of unknown displacement gradients in the element, which must be determined as part of the solution.

4.      Similarly, the virtual displacement gradient is written as

$\frac{\partial \delta v_i}{\partial x_j}=\left({\displaystyle \sum _{a=1}^{N_e}\frac{\partial N^a(\xi )}{\partial {\xi }_m}}\delta v_i^a+{\displaystyle \sum _{k=1}^p\frac{J(0)}{J(\xi )}}\delta {\alpha }_i^{(k)}{\delta }_{km}{\xi }_k\right)\frac{\partial {\xi }_m}{\partial x_j}$

where $\delta {\alpha }_i^{(k)}$ is a variation in the internal displacement gradient field for the element.  

5.      These expressions are then substituted into the virtual work equation, which must now be satisfied for all possible values of virtual nodal displacements $\delta v_i^a$ and virtual displacement gradients $\delta {\alpha }_i^{(k)}$.  At first sight this procedure appears to greatly increase the size of the global stiffness matrix, because a set of unknown displacement gradient components must be calculated for each element.  However, the unknown ${\alpha }_i^{(k)}$ are local to each element, and can be eliminated while computing the element stiffness matrix.  The procedure to do this can be shown most clearly in a sample code.

 

A sample small-strain, linear elastic code with incompatible mode elements is provided in Femlinelast_incompatible_modes.mws.  When run with the input file shear_locking_demo.txt it produces the results shown in the figure.  The incompatible modes clearly give a spectacular improvement in the performance of the element.

 

HEALTH WARNINGS: (i) The procedure outlined here only works for small-strain problems.  Finite strain versions exist but are somewhat more complicated.  (ii)  Adding strain variables to elements can dramatically improve their performance, but this procedure must be used with great care to ensure that the strain and displacement degrees of freedom are independent variables. For further details see Simo,  J. C., and M. S. Rifai, Int J. Num. Meth in Eng.  29, pp. 1595$–$1638, 1990. Simo,  J. C., and F. Armero, Int J. Num. Meth in Eng., 33, pp. 1413$–$1449, 1992.

 

 

8.6.2 Volumetric locking and reduced integration elements

 

Volumetric locking can be illustrated using a simple boundary value problem. Consider a long hollow cylinder with internal radius a and external radius b as shown in the figure.  The solid is made from a linear elastic material with Young’s modulus $E$ and Poisson’s ratio $\nu$.  The cylinder is loaded by an internal pressure $p_a$ and deforms in plane strain. 

The analytical solution to this problem is given in Section 4.1.9.

 

The figures below compare the analytical solution to a finite element solution with standard 4 noded plane strain quadrilateral elements.  Results are shown for two values of Poisson’s ratio $\nu$.  The dashed lines show the analytical solution, while the solid line shows the FEA solution.

 

           

The two solutions agree well for $\nu =0.3$, but the finite element solution grossly underestimates the displacements as Poisson’s ratio is increased towards 0.5 (recall that the material is incompressible in the limit $\nu =0.5$ ).  In this limit, the finite element displacements tend to zero $–$ this is known as `volumetric locking’

 

The error in the finite element solution occurs because the finite element interpolation functions are unable to properly approximate a volume preserving strain field. In the incompressible limit, a nonzero volumetric strain at any of the integration points gives rise to a very large contribution to the virtual power.  The interpolation functions can make the volumetric strain vanish at some, but not all, the integration points in the element. 

 

Volumetric locking is a much more serious problem than shear locking, because it cannot be avoided by refining the mesh.  In addition, all the standard fully integrated finite elements will lock in the incompressible limit; and some elements show very poor performance even for Poisson’s ratios as small as $0.45$.  Fortunately, most materials have Poisson’s ratios around 0.3 or less, so the standard elements can be used for most linear elasticity and small-strain plasticity problems.   To model rubbers, or to solve problems involving large plastic strains, the elements must be redesigned to avoid locking.

 

Reduced Integration is the simplest way to avoid locking.   The basic idea is simple: since the fully integrated elements cannot make the strain field volume preserving at all the integration points, it is tempting to reduce the number of integration points so that the constraint can be met.  `Reduced integration’ usually means that the element stiffness is integrated using an integration scheme that is one order less accurate than the standard scheme.   The number of reduced integration points for various element types is listed in the table below.  The coordinates of the integration points are listed in the tables in Section 8.1.12.

 

 

HEALTH WARNING: Notice that the integration order cannot be reduced for the linear triangular and tetrahedral elements $–$ these elements should not be used to model near incompressible materials, although in desperation you can a few such elements in regions where the solid cannot be meshed using quadrilaterals.

 

Remarkably, reduced integration completely resolves locking in some elements (the quadratic quadrilateral and brick), and even improves the accuracy of the element.  As an example, the figure below shows the solution to the pressurized cylinder problem, using both full and reduced integration for 8 noded quadrilaterals.  With reduced integration, the analytical and finite element results are indistinguishable.

                

 

 

Reduced integration does not work in 4 noded quadrilateral elements or 8 noded brick elements.  For example, the figure on the right shows the solution to the pressure vessel problem with linear (4 noded) quadrilateral elements with reduced integration.  The solution is clearly a disaster.  The error occurs because the stiffness matrix is nearly singular $–$  the system of equations includes a weakly constrained deformation mode. This   phenomenon is known as `hourglassing’ because of the characteristic shape of the spurious deformation mode. 

 

Selectively Reduced Integration can be used to cure hourglassing.   The procedure is illustrated most clearly by modifying the formulation for static linear elasticity. To implement the method: 

1.      The volume integral in the virtual work principle is separated into a deviatoric and volumetric part by writing

${\displaystyle \underset{R}{\int }{\sigma }_{ij}\delta {\epsilon }_{ij}d{V}_0}={\displaystyle \underset{R}{\int }\left({\sigma }_{ij}\delta {\epsilon }_{ij}-\frac{{\sigma }_{kk}}{3}\delta {\epsilon }_{qq}\right)d{V}_0}+{\displaystyle \underset{R}{\int }\frac{{\sigma }_{kk}}{3}\delta {\epsilon }_{qq}d{V}_0}$

 Here, the first integral on the right hand side vanishes for a hydrostatic stress.

2.      Substituting the linear elastic constitutive equation and the finite element interpolation functions into the virtual work principle, we find that the element stiffness matrix can be reduced to

$k_{aibk}^{(l)}={\displaystyle \underset{V_e^{(l)}}{\int }\left(C_{ijkl}\frac{\partial N^a(x)}{\partial x_j}\frac{\partial N^b(x)}{\partial x_l}-\frac{1}{3}{C}_{ppkl}\frac{\partial N^a(x)}{\partial x_i}\frac{\partial N^b(x)}{\partial x_l}\right)dV}+{\displaystyle \underset{V_e^{(l)}}{\int }\frac{1}{3}{C}_{ppkl}\frac{\partial N^a(x)}{\partial x_i}\frac{\partial N^b(x)}{\partial x_l}dV}$

3.      When selectively reduced integration is used, the first volume integral is evaluated using the full integration scheme; the second integral is evaluated using reduced integration points.

 

Selective reduced integration has been implemented in the sample program fem_selective_reduced_integration.mws.  When this code is run with the input file volumetric_locking_demo.txt it produces the results shown in the figure.  The analytical and finite element solutions agree, and there are no signs of hourglassing.

 

In many commercial codes, the `fully integrated’ elements actually use selective reduced integration.

 

 

The ‘B-bar’ method:  Like selective reduced integration, the B-bar method works by treating the volumetric and deviatoric parts of the stiffness matrix separately.  Instead of separating the volume integral into two parts, however, the B-bar method modifies the definition of the strain in the element.  Its main advantage is that the concept can easily be generalized to finite strain problems.   Here, we will illustrate the method by applying it to small-strain linear elasticity.  The procedure starts with the usual virtual work principle 

${\displaystyle \underset{R}{\int }{\sigma }_{ij}[{\epsilon }_{ij}(u_k)]\delta {\epsilon }_{ij}dV-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}dV-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}\delta v_{i}dA=0}}}$

In the B-bar method

1.      We introduce a new variable to characterize the volumetric strain in the elements.  Define

$\omega =\frac{1}{3V_e}{\displaystyle \underset{V_e}{\int }{\epsilon }_{kk}dV}=B_{bk}^{(vol)}{u}_k^b{B}_{bk}^{(vol)}=\frac{1}{3V_e}{\displaystyle \underset{V_e}{\int }\frac{\partial N^b}{\partial x_k}dV}$

where the integral is taken over the volume of the element.

2.      The strain variation in each element is replaced by the approximation ${\overline{\epsilon }}_{ij}={\epsilon }_{ij}+(\omega -{\epsilon }_{kk}/3){\delta }_{ij}$

3.      Similarly, the virtual strain in each element is replaced by $\delta {\overline{\epsilon }}_{ij}=\delta {\epsilon }_{ij}+(\delta \omega -\delta {\epsilon }_{kk}/3){\delta }_{ij}$

This means that the volumetric strain in the element is everywhere equal to its mean value.  The virtual work principle is then written in terms of ${\overline{\epsilon }}_{ij}$ and $\delta {\overline{\epsilon }}_{ij}$ as

${\displaystyle \underset{R}{\int }{\sigma }_{ij}[{\overline{\epsilon }}_{ij}(u_k)]\delta {\overline{\epsilon }}_{ij}dV-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}dV-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i^{*}\delta v_{i}dA=0}}}$

Finally, introducing the finite element interpolation and using the constitutive equation yields the usual system of linear finite element equations, but with a modified element stiffness matrix given by

$k_{aibk}^{(l)}={\displaystyle \underset{V_e^{(l)}}{\int }{C}_{pjql}{\overline{B}}_{pji}^a{\overline{B}}_{qlk}^{b}dV}{\overline{B}}_{pji}^a=\frac{\partial N^a}{\partial x_j}{\delta }_{ip}+\left(B_{ia}^{(vol)}-\frac{1}{3}\frac{\partial N^a}{\partial x_i}\right){\delta }_{pj}$

This expression can be integrated using a standard, full integration scheme.

 

The B-bar method has been implemented in the sample code FEM_Bbar.mws.  The code can be run with the input file volumetric_locking_demo.txt.  Run the code for yourself to verify that the analytical and finite element solutions agree, and there are no signs of hourglassing.

Hourglass base vectors
Linear quadrilateral	${\Gamma }^{a(1)}=(+1,-1,+1,-1)$
Linear brick	$\begin{array}{l}{\Gamma }^{a(1)}=(+1,+1,-1,-1,-1,-1,+1,+1)\\ {\Gamma }^{a(2)}=(+1,-1,-1,+1,-1,+1,+1,-1)\\ {\Gamma }^{a(3)}=(+1,-1,+1,-1,+1,-1,+1,-1)\\ {\Gamma }^{a(4)}=(-1,+1,-1,+1,+1,-1,+1,-1)\end{array}$

 

Reduced integration with hourglass control: Hourglassing in 4 noded quadrilateral and 8 noded brick elements can also be cured by adding an artificial stiffness to the element that acts to constrain the hourglass mode.   The stiffness must be carefully chosen so as to influence only the hourglass mode.  Only the final result will be given here $–$ for details see D.P. Flanagan and T. Belytschko, International J. Num Meth in Engineering, 17, pp. 679$–$706, (1981).  To compute the corrective term:

1.      Define a series of `hourglass base vectors’  ${\Gamma }^{a(i)}$ which specify the displacements of the ath node in the ith hourglass mode.   The 4 noded quadrilateral element has only one hourglass mode; the 8 noded brick has 4 modes, listed in the table.

2.      Calculate the `hourglass shape vectors’ for each mode as follows

${\gamma }^{a(i)}={\Gamma }^{a(i)}-\frac{\partial N^a(\xi =0)}{\partial x_j}{\displaystyle \sum _{b=1}^{n_e}{\Gamma }^{b(i)}{x}_j^b}$

where $n_e$ denotes the number of nodes on the element.

3.      The modified stiffness matrix for the element is written as

$k_{aibk}^{(l)}={\displaystyle \underset{V_e^{(l)}}{\int }{C}_{ijkl}\frac{\partial N^a(x)}{\partial x_j}\frac{\partial N^b(x)}{\partial x_l}dV}+\kappa V_e^{(l)}{\displaystyle \sum _m{\gamma }^{a(m)}{\gamma }^{b(m)}}$

where $V_e$ denotes the volume of the element, and $\kappa$ is a numerical parameter that controls the stiffness of the hourglass resistance.  Taking $\kappa =0.01\mu \left(\partial N^a/\partial x_j\right)\left(\partial N^a/\partial x_j\right)$ where $\mu$ is the elastic shear modulus works well for most applications.  If $\kappa$ is too large, it will seriously over-stiffen the solid.

 

Sample code: Reduced integration with hourglass control has been implemented in the sample code Fem_hourglasscontrol.mws.   When run with the input file volumetric_locking_demo.txt it produces the results shown in the picture.  Hourglassing has clearly been satisfactorily eliminated.

 

HEALTH WARNING: Hourglass control is not completely effective: it can fail for finite strain problems and can also cause problems in a dynamic analysis, where the low stiffness of the hourglass modes can introduce spurious low frequency vibration modes and low wave speeds.

 

 

8.6.3 Hybrid elements for modeling near-incompressible materials

 

The bulk elastic modulus is infinite for a fully incompressible material, which leads to an infinite stiffness matrix in the standard finite element formulation (even if reduced integration is used to avoid locking).   This behavior can cause the stiffness matrix for a nearly incompressible material to become ill-conditioned, so that small rounding errors during the computation result in large errors in the solution.

 

Hybrid elements are designed to avoid this problem.  They work by including the hydrostatic stress distribution as an additional unknown variable, which must be computed at the same time as the displacement field.  This allows the stiff terms to be removed from the system of finite element equations.  The procedure is illustrated most easily using isotropic linear elasticity, but in practice hybrid elements are usually used to simulate rubbers or metals subjected to large plastic strains.

 

Hybrid elements are based on a modified version of the principle of virtual work, as follows. The virtual work equation (for small strains) is re-written as

${\displaystyle \underset{R}{\int }\left(S_{ij}\delta {\epsilon }_{ij}\right)d{V}_0}+{\displaystyle \underset{R}{\int }p\delta {\epsilon }_{qq}d{V}_0}+{\displaystyle \underset{R}{\int }\delta p\left(\frac{{\sigma }_{kk}}{3K}-\frac{p}{K}\right)d{V}_0}-{\displaystyle \underset{R}{\int }{b}_i\delta v_{i}d{V}_0-{\displaystyle \underset{{\partial }_{2}R}{\int }{t}_i\delta v_{i}dA}}=0$

Here

1.      $S_{ij}={\sigma }_{ij}-{\sigma }_{kk}{\delta }_{ij}/3$ is the deviatoric stress, determined from the displacement field;

2.      ${\sigma }_{kk}/3$ is the hydrostatic stress, again determined from the displacement field;

3.      p is an additional degree of freedom that represents the (unknown) hydrostatic stress in the solid;

4.      $\delta p$ is an arbitrary variation in the hydrostatic stress;

5.      $K=(1/3)\partial {\sigma }_{kk}/\partial {\epsilon }_{pp}$ is the bulk modulus of the solid.

The modified virtual work principle states that, if the virtual work equation is satisfied for all kinematically admissible variations in displacement $\delta v_i$ and strain $\delta {\epsilon }_{ij}=\left(\partial \delta v_i/\partial x_j+\partial \delta v_j/\partial x_i\right)$ and all possible variations in pressure $\delta p$, the stress field will satisfy the equilibrium equations and traction boundary conditions, and the pressure variable p will be equal to the hydrostatic stress in the solid.

  

The finite element equations are derived in the usual way

1.      The displacement field, virtual displacement field and position in the each element are interpolated using the standard interpolation functions defined in Sections 8.1.9 and 8.1.10  as

$u_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)u_i^a}\delta v_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)\delta v_i^a}{x}_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)x_i^a}$

2.      Since the pressure is now an independent variable, it must also be interpolated.  We write

$p(x)={\displaystyle \sum _{a=1}^{N_p}{M}^a({\xi }_j)p^a}\delta p(x)={\displaystyle \sum _{a=1}^{N_p}{M}^a({\xi }_j)\delta p^a}{x}_i={\displaystyle \sum _{a=1}^{N_e}{N}^a({\xi }_j)x_i^a}$

where $p^a$ are a discrete set of pressure variables, $\delta p^a$ is an arbitrary change in these pressure variables, $M^a$ are a set of interpolation functions for the pressure, and $N_p$ is the number of pressure variables associated with the element.  The pressure need not be continuous across neighboring elements, so that independent pressure variables can be added to each element.  The following schemes are usually used

a.       In linear elements (the 3 noded triangle, 4 noded quadrilateral, 5 noded tetrahedron or 8 noded brick) the pressure is constant.  The pressure is defined by its value at the centroid of each element, and the interpolation functions are constant.

b.      In quadratic elements (6 noded triangle, 8 noded quadrilateral, 10 noded tetrahedron or 20 noded brick) the pressure varies linearly in the element.  Its value can be defined by the pressure at the corners of each element, and interpolated using the standard linear interpolation functions.

3.      For an isotropic, linear elastic solid with shear modulus $\mu$ and Poisson ratio $\nu$ the deviatoric stress is related to the displacement field by  $S_{ij}=\mu ({\delta }_{il}{\delta }_{jk}+{\delta }_{ik}{\delta }_{jl}-2{\delta }_{ij}{\delta }_{kl}/3)\left(\partial u_k/\partial x_l\right)$, while the hydrostatic stress is ${\sigma }_{kk}/3=K(\partial u_k/\partial x_k)$, where $K=E/3(1-2\nu )$ is the bulk modulus.

4.      Substituting the linear elastic equations and the finite element interpolation functions into the virtual work principle leads to a system of equations for the unknown displacements $u_k^b$ and pressures $p^b$ of the form

$\begin{array}{l}{K}_{aibk}{u}_k^b+Q_{aib}{p}^b=F_i^a\forall \{a,i\}:x_k^a\text{not on }{\partial }_{\text{1}}R\\ Q_{akb}{u}_k^b-{\Pi }_{ab}{p}^b=0\\ u_i^a=u_i^{*}(x_i^a)\text{}\forall \{a,i\}:x_k^a\text{on }{\partial }_{\text{1}}R\end{array}$

where the global stiffness matrices $K,Q,\Pi$ are obtained by summing the following element stiffness matrices

$\begin{array}{l}{k}_{aibk}={\displaystyle \underset{V_e}{\int }{C}_{ijkl}\frac{\partial N^a(x)}{\partial x_j}\frac{\partial N^b(x)}{\partial x_l}-K\frac{\partial N^a(x)}{\partial x_i}\frac{\partial N^b(x)}{\partial x_k}dV}{q}_{aib}={\displaystyle \underset{V}{\int }\frac{\partial N^a(x)}{\partial x_i}{M}^b(x)dV}\\ {\pi }_{ab}={\displaystyle \underset{V_e}{\int }\frac{1}{K}{M}^a(x)M^b(x)dV}\end{array}$

and the force F is defined in the usual way.  The integrals defining $k_{iakb}$ may be evaluated using the full integration scheme or reduced integration (hourglass control may be required in this case).  The remaining integrals must be evaluated using reduced integration to avoid element locking.

5.      Note that, although the pressure variables are local to the elements, they cannot be eliminated from the element stiffness matrix, since doing so would reduce the element stiffness matrix to the usual, non-hybrid form.  Consequently, hybrid elements increase the cost of storing and solving the system of equations.