Mixed displacement–pressure formulations and suitable finite elements for multimaterial problems with compressible and incompressible models

Multimaterial problems in linear and nonlinear elasticity are some of the least explored using mixed finite element formulations with higher-order elements. The fundamental issue in adapting the mixed displacement–pressure formulations with linear and higher-order continuous elements for the pressure field is their inability to capture pressure and stress jumps across material interfaces. In this paper, for the first time in literature, we perform comprehensive studies of multimaterial problems in elasticity consisting of compressible and incompressible material models using the mixed displacement–pressure formulation to assess the performance of different element types in accurately resolving pressure fields within the domains and pressure jumps across material interfaces. In particular, inf–sup stable displacement–pressure combinations with element-wise discontinuous pressure for triangular and tetrahedral elements are considered and their performance is assessed along with the Q1/P0 element and Taylor– Hood elements using several numerical examples. The results show that Taylor–Hood elements fail to capture the stress jumps due to the continuity of DOFs across elements, the Crouzeix– Raviart (P2b/P1dc) element yields substantially poor pressure fields despite a significant increase in pressure degrees of freedom and that the P3/P1dc element produces superior quality results fields when compared with the P2b/P1dc element.


Introduction
Elasticity problems consisting of multiple materials are ubiquitous in nature, e.g., biological tissues and man-made composites.For rubber-like polymers and soft tissues, the stress-strain relation is defined by hyperelastic models, which also form the building blocks for the constitutive relations for soft multifunctional composites such as electroactive polymers [1], magnetoactive polymers [2], hydrogels [3], etc.The typical deformation behaviour of materials, in general, can range from the compressible regime with significant volumetric changes to the incompressible regime with zero volumetric change.
While compressible materials are relatively straightforward to simulate using pure displacement formulations, incompressible hyperelastic models require sophisticated finite element formulations for computing accurate numerical results.Formulations such as  -bar formulation [4], average nodal strain formulation [5], reduced integration method with hourglass control [6], enhancedstrain method [7], enhanced assumed strain (EAS) methods [8], energy-sampling stabilisation [9], smoothed FEM [10][11][12] etc., that take the contribution from the volumetric part of the strain energy function into the stiffness matrix end up being computationally expensive due to the increased ill-conditioning of the resulting stiffness matrix as the incompressibility limit is approached, that is, as Poisson's ratio approaches 0.5.For simulating the truly incompressible material models accurately, mixed displacement-pressure formulation, also referred to as the hybrid formulation, is the only viable option.
The majority of mixed displacement-pressure formulations with inf-sup stable elements such as Taylor-Hood (P2/P1 and Q2/Q1) elements [13,38] and Bézier family [16], subdivision stabilised NURBS elements [15], and mixed stabilised formulations [19,24,27], have been studied extensively for problems with a single material.However, problems of practical interest often consist of multiple materials, requiring numerical methods that are applicable to such problems.A fundamental requirement in simulating solid mechanics problems with multiple materials is accurately capturing the stress discontinuities across material interfaces.While the mixed displacement-pressure formulation with the Q1/P0 element (linear quadrilateral/hexahedral for displacement and elementwise constant for pressure) is readily applicable to multimaterial problems, higher-order elements are not so readily adaptable for problems with multiple materials.Although the Taylor-Hood family of elements satisfy the inf-sup stability condition, their straightforward applicability is limited to incompressible solid mechanics problems with a single material.For multi-material problems, direct application of the Taylor-Hood elements results in incorrect pressure (and hence stress) values at the material interfaces.This is due to the continuity of pressure degrees of freedom (DOFs) across elements, which makes it impossible to capture pressure (and hence stress) jumps at material interfaces without additional effort.
To accurately resolve pressure jumps at the material interfaces using linear or higher-order elements for pressure, the pressure DOFs for the elements sharing a material interface must be duplicated, as done in Kadapa et al. [3].However, such an approach of duplicating pressure DOFs is only possible to implement with knowing the additional information, such as the area/volume the element is attached to.Note that we cannot prepare the duplicate nodes based on the material information alone, as disjointed locations can share the same material number, for example, particles distributed in a matrix.Moreover, the approach of duplication of pressure DOFs is challenging to implement in parallel codes for distributed memory architectures based on a message-passing interface; this requires a significant amount of data transfer among processes.An alternative is to generate non-matching meshes for each component, but this requires a contact formulation to enforce the displacement continuity, which adds unnecessary overhead.Therefore, pressure spaces must be discontinuous to accurately resolve the pressure jumps across material interfaces.
Although the well-known Q1/P0 element [17,39] is readily adaptable to such problems, its lack of inf-sup stability and lowerorder nature forces one to use extremely finer meshes for accurate results [1].Among higher-order elements with discontinuous pressure spaces, the Crouzeix-Raviart element (P2b/P1dc -quadratic triangular/tetrahedral element with a cubic bubble for displacements and linear element-wise discontinuous space for pressure) [40][41][42] and P3/P1dc element (cubic triangular/tetrahedral element for displacements and linear element-wise discontinuous space for pressure) [42,43] are of particular interest because they allow for the easy of generation of triangular and tetrahedral meshes for complex geometries.However, these elements have not been explored in the context of nonlinear elasticity problems.Motivated by the need to accurately simulate nonlinear elasticity problems consisting of multiple materials of which one or all materials are incompressible, in this work, we evaluate the accuracy and performance of P2b/P1dc and P3/P1dc elements in the context of linear and nonlinear elasticity problems, particularly those made up multiple materials, for the first time in literature to the best of our knowledge.The accuracy of results is evaluated by calculating error norms for the examples of thick-walled cylinders and spheres made of linear elastic materials and subjected to internal pressure.The performance of the elements is further demonstrated by studying benchmark examples and their extensions in nonlinear elasticity with hyperelastic materials.

Kinematics, strain and stress measures and constitutive relations
We consider an arbitrary solid body with  0 as the reference configuration and   as the current configuration the body assumes under the action of external loads.A nonlinear deformation map  ∶  0 →   takes the point  in the reference configuration,  0 , to the point  in the current configuration,   , so that the displacement is  =  − .Assuming that the map  is smooth and invertible so that the deformation gradient  is well defined. and its determinant  are given as where  is the second-order identity tensor.The right Cauchy-Green deformation tensor, , the Green-Lagrange strain tensor, , and the infinitesimal strain tensor, , are defined as In the case of truly incompressible materials in the finite strain regime, the deformation of the solid is such that the total volume change at any point in the domain is zero.This acts as a constraint on the displacement field that needs to be solved in the elasticity problem.Mathematically, the incompressibility constraint in the finite strain regime is given as,  = 1. (5) To model hyperelastic materials in the incompressible finite strain regime, the deformation gradient is decomposed into deviatoric and volumetric components as where  vol ∶=  1∕3 , and  dev ∶=  −1∕3  , (7) and the modified version of the deformation gradient is defined as Accordingly, the modified versions of the right Cauchy-Green tensor and the Green-Lagrange strain tensors are given as The nonlinear behaviour of soft polymers and soft biological tissues is defined by using hyperelastic models in the form of strain energy density functions.For nearly incompressible materials, the strain energy functions for hyperelastic models are assumed to be split into deviatoric and volumetric parts as, A wide variety of functions for deviatoric and volumetric are available in the literature.In this work, Neo-Hookean and Mooney-Rivlin models are considered for the deviatoric part.

Mixed displacement-pressure formulation in small strains
In this work, the mixed formulation in small strains is tuned for linear isotropic elastic materials.For the mixed displacementpressure formulation in small strains, the infinitesimal strain tensor is decomposed into deviatoric and volumetric parts as with ,  (14) where the subscripts ,  and  denote the coordinate components, and the repeated index means summation.For example, .With  as the additional unknown, the Cauchy stress for linear elastic material for the mixed formulation is calculated as (15) and in component form as For compressible materials ( < 0.5), hydrostatic pressure, , is related to the volumetric strain via the relation where  is the bulk modulus.For incompressible materials,  acts as a Lagrange multiplier enforcing the volumetric constraint Considering both compressible and incompressible materials, the total potential energy for the mixed displacement-pressure formulation in small strains can be written as where is the energy contribution from the external body force,  0 , and traction,  0 .Following the variational calculus, the first variation (()) and second variation (()) of  are where is the fourth-order elasticity tensor.
Using the finite element approximations for the solution variables and their variations, where   and   are the basis functions, respectively, for the displacement and pressure fields, and  and , respectively, are the nodal degrees of freedom (DOFs), respectively, for the displacement and pressure fields, the coupled matrix system for the mixed formulation can be written as where where  is the fourth-order elasticity tensor (C) expressed in the matrix form. 0 and  0 are the gradient-displacement and divergence-displacement matrices with respect to the reference configuration, which, for a single basis function,   , are given as

Mixed displacement-pressure formulation in finite strains
We adapt the mixed displacement-pressure formulation that is recently proposed by Kadapa and Hossain [17] because of its advantage in yielding symmetric global stiffness matrices irrespective of the volumetric energy function.Following Kadapa and Hossain [17], the total energy functional for the mixed displacement-pressure formulation in the finite strain regime is given as, where   is the generic energy functional that accounts for the effect of volumetric deformation in the compressible as well as incompressible regime while considering arbitrary volumetric energy function.  enforces the incompressibility constraint  =1 using the Lagrange multiplier approach for the truly incompressible case while it uses a more generic constraint by incorporating the volumetric energy function  vol [17].For the generic case, we have   ( , ) in the perturbed Lagrangian form as Computer Methods in Applied Mechanics and Engineering 432 (2024) 117354 with  as the Lagrange multiplier.Ĵ and θ are calculated using with  • = det ( ( • )).The subscript  in Eq. ( 32) denotes the previously converged load step.From the generic case in (31), we can recover the truly incompressible case by setting Ĵ = 1 and θ = 0 based on the user input.The reader is referred to Kadapa and Hossain [17] for the comprehensive details on the mixed formulation and its demonstrated advantages using benchmark numerical examples.
By taking the first variation () of the total energy functional given in Eq. ( 30), we get Using the finite element approximations for the displacement and pressure fields and their variations in (24), the semi-discrete form for the mixed formulation in finite strains can be written as, where  is the total stress and discrete gradient-displacement matrix with respect to the current configuration () for a single basis function is given as We use the Newton-Raphson scheme to solve the set of nonlinear Eqs.(34a) and (34b) in an efficient manner.With  + 1 and  denoting the current and previously converged load steps, respectively, and  + 1 and  denoting the current and previous iterations at the current load step, the displacement DOFs,  +1 , and pressure DOFs,  +1 , at the current iteration of the current load step are computed as the iterative increments,  and , from the respective quantities at the previous iteration,  ()  +1 and  () +1 , as By adapting the Newton-Raphson scheme, we arrive at the coupled matrix system given as where where  is material tangent tensor (e) expressed in the matrix form and  is the divergence-displacement matrix with respect to the current configuration, which, for a single basis function,   , is given as Note that   contains contributions from geometric and material nonlinearities.Following the derivation above, the effective first Piola-Kirchhoff stress tensor ( ) for the mixed formulation is defined as Computer Methods in Applied Mechanics and Engineering 432 (2024) 117354 Fig. 1.A typical interface between two domains made of materials with different material properties.
In the above equation, deviatoric component of the first Piola-Kirchhoff stress tensor,  is computed solely from the deviatoric part of the strain energy function,  dev .From  , the Cauchy stress tensor becomes, The material tangent tensor (e) of order four is computed as

Finite element spaces
It is well-established that the finite element spaces for the displacement and pressure for the mixed formulation cannot be chosen arbitrarily [41].The combination of displacement-pressure spaces must satisfy the inf-sup condition, and a pressure space that is not inf-sup compatible with the displacement space yields spurious oscillations in the pressure field.Because of this restriction, the number of such spaces for the displacement-pressure combination is limited.Among the class of finite elements that satisfy the inf-sup condition, Taylor-Hood elements such as P2/P1, Q2/Q1, and their higher-order extensions are popular.Recently, their extensions to NURBS [15] and Bézier elements, e.g., BT2/BT1, BQ2/BQ1 and BW2/BW1 [16,44] have also been proposed.Among the elements that do not satisfy the inf-sup condition, the Q1/P0 element is widely used because of its simplicity despite yielding a checkerboard pattern for pressure.
However, Taylor-Hood elements pose some serious issues when applied to problems with material interfaces.To understand this, let us consider a material interface between two domains with different material properties as shown in Fig. 1.The Young's moduli of the two domains are  1 and  2 , and their respective Poisson's ratio are  1 and  2 .The displacement and pressure fields for the first domain are  1 and  1 , respectively, and the corresponding quantities for the second domain are  2 and  2 .For the mixed displacement-pressure formulation, we can express the Cauchy stresses in the th domain as The displacement field should be continuous at the interface for the perfectly bonded interface.Moreover, the Euler and Cauchy stress principle [45] mandates the continuity of tractions at the interface.The continuity of displacements is easily ensured by the  0 continuity of the adapted Lagrange family of elements for the displacement field.But, continuity of tractions ( () = − (−) ), as given by requires that the pressure fields  1 and  2 must be discontinuous across the interface.However, because of the  0 continuity of pressure DOFs for the Taylor-Hood elements, to capture such discontinuities in the pressure field either the pressure DOFs at the interface must be duplicated [3], which is quite cumbersome for generic problems in 3D, or a suitable contact formulation must be employed if both displacement and pressure DOFs are duplicated or non-matching meshes are used, which adds unnecessary overhead.As an alternative, we propose elements with discontinuous pressure spaces.
In this work, we consider the Q1/P0 element, the P2/P1 element, the Crouzeix-Raviart (P2b/P1dc) element [40], and the P3/P1dc element.The P2b/P1dc element is a quadratic triangular/tetrahedral element enriched with a cubic bubble for the displacement (P2b) and element-wise linear discontinuous element (P1dc) for pressure.In the P3/P1dc element, the displacement field is approximated with cubic triangular/tetrahedral element and the pressure is approximated with element-wise linear discontinuous spaces.The displacement and pressure nodes for the different 2D elements considered in this work are illustrated in Fig. 2.

Numerical examples
We study the accuracy of Q1/P0, P2/P1, P2b/P1dc and P3/P1dc elements using several benchmark examples especially in nonlinear elasticity and compare their relative merits and issues.The formulations are implemented in an in-house code written in C++.The meshes are generated using GMSH [46] and the results are postprocessed using ParaView [47].
The strain energy functions considered in this work are given below.
• Neo-Hookean model • Mooney-Rivlin model • Volumetric energy function where  is the shear modulus,  is the bulk modulus,  01 and  10 are the material constants.

Thick-walled cylinder under internal pressure
In the first example, we consider a thick-walled cylinder subjected to internal pressure.The 2D model is assumed to be in plane-strain condition.In addition to the simple thick-walled cylinder made of a single material, see Fig. 3(a), we also consider a composite cylinder that is made of two different materials, see Fig. Since the condensed form of Q1/P0 is not applicable to perfectly incompressible cases, Poisson's ratio for the Q1/P0 element is taken as 0.4999.The analytical solutions are provided in Appendix.

Simple thick-walled cylinder
Convergence analysis is conducted with five uniform meshes shown in Fig. 4. Graphs of  2 relative error norms in displacement, pressure and stress against mesh refinement are shown in Figs. 5 and 6, respectively, for  = 0.3 and  = 0.5.As shown, all four element types considered in the study yield the expected convergence rates.A closer look at the convergence rates and the magnitude  of errors reveals a surprising outcome in the convergence of the pressure field.While the magnitude of errors in the displacement of the P2b/P1dc element is only slightly lower than that of the P2/P1 element, the magnitude of error in the pressure field obtained with the P2b/P1dc element is even higher than that of the Q1/P0 element, and at least an order of magnitude higher than that of the P2/P1 element.The error magnitude of the pressure field obtained with the P2b/P1dc element is almost two orders higher and the error magnitude of stress is one order of magnitude higher when compared with the P3/P1dc element, indicating a subpar performance of the P2b/P1dc element despite a substantial increase in the number of pressure DOFs when compared with the Q1/P0 and P2/P1 elements and while having the same number of pressure DOFs as that of the P3/P1dc element.The performance of the P3/P1dc element is substantially superior to that of the P2b/P1dc element although both elements have the same number of pressure DOFs.
The subpar performance of the P2b/P1dc element can be further understood through warped plots of the pressure field obtained.As shown in Fig. 7, the pressure field obtained with the P2b/P1dc element is erratic as opposed to the smooth pressure fields obtained with the P2/P1 and P3/P1dc elements.

Composite thick-walled cylinder
Similar convergence studies are conducted for the thick-walled cylinder made of two materials, referred to as the composite cylinder henceforth, using five meshes shown in Fig. 8. Graphs of error norms are shown in Figs. 9 and 10 for two different values of Poisson's ratio for the same values of Young's moduli,  1 = 200 MPa and  2 = 20 MPa.The graphs show the expected rate of convergence for displacement, pressure and stress fields for the Q1/P0, P2b/P1dc and P3/P1dc elements.However, the rate of convergence for the P2/P1 element deteriorates substantially in all the fields: the displacement convergence rate is one instead of three and pressure and stress converge at a rate of 0.5 instead of two.This drop in convergence rates for the P2/P1 element is due to the fact that the pressure nodes (and hence DOFs) are continuous across the elements.Therefore, the default implementation of the P2/P1 element is incapable of capturing the pressure jumps across material interfaces, see Fig. 11(b).On the contrary, Q1/P0, P2b/P1dc and P3/P1dc elements capture the pressure jumps accurately.However, within each cylindrical domain, the pressure field obtained with the P2b/P1dc element (see Fig. 11(c)) shows the same erratic behaviour as observed in the case of the simple     cylinder.The pressure field obtained with the P3/P1dc element is not only discontinuous at the material interface but also smoother within each domain, as shown in Fig. 11(d).
At this point, it is worth highlighting the superior accuracy of the P3/P1dc element in all the fields.With only three additional nodes for the displacement field when compared with the P2b/P1dc element, the convergence rates of displacement and pressure (and stress) increase, respectively, by an order of one and 0.5, and at the same time the magnitudes of errors in all the quantities are lower by at least one order when compared with the P2b/P1dc element.The additional computational cost of the P3/P1dc element is worth it when we consider the accuracy of the results obtained, particularly pressure and stress fields.

Thick-walled composite sphere under internal pressure
We consider the thick-walled composite sphere to study the performance 3D versions of P2/P1, P2b/P1dc and P3/P1dc elements.The radii of the spheres are 100, 150 and 200 mm.Due to the symmetry, only 1/8th of the sphere with symmetry boundary conditions is considered for the analysis.The Young's moduli of the inner and outer cylinders, respectively, are  1 = 200 MPa and  2 = 20 MPa, and the Poisson's ratio is  1 =  2 = 0.5.Four meshes shown in Fig. 12 are used for the analysis, and the graphs of displacement, pressure and stress fields are shown in Fig. 13.As shown, the error norms for all three element types show the   same trend as observed in the 2D case.The P2/P1 element fails to capture the pressure jumps accurately, and the P3/P1dc element shows superior performance when compared with the P2b/P1dc element, yielding a higher convergence rate and lower magnitude of errors.The contour plots of pressure presented in Fig. 14 for mesh 2 show the inability of the P2/P1 element to accurately capture pressure jumps at the interface, and smoother pressure fields along with accurate pressure jumps obtained with the P3/P1dc element.

Cook's membrane in 2D
We now consider the original and modified cases of the popular Cook's membrane benchmark.The first case is the Cook's membrane problem with a single material, and the second case is the modified one with two materials with different material properties, as shown in Fig. 15.

Cook's membrane with single material
For this case, Young's modulus is  = 240.565MPa and Poisson's ratio is  = 0.4999.As shown in Fig. 16, five meshes are used for the analysis using triangular elements.Mapped meshes with corresponding number of elements per side are used for the Q1/P0   element.Only five load steps are used for the simulations with the Neo-Hookean model.The Y-displacement of point A with mesh refinement using different element types is shown in Fig. 17 for linear elastic and Neo-Hookean materials.All the element types converge with mesh refinement and match well with the reference values from the literature.Although the solution with the Q1/P0 element converges with mesh refinement, the pressure field obtained shows the well-known checkerboard pattern, see Fig. 18.The pressure field obtained with P2/P1 and P3/P1dc elements is smoother than that obtained with the P2b/P1 element, consistent with the results obtained in the thick-walled cylinder example.

Cook's membrane with two materials
For the two-material version, two cases are considered.The Young's modulus and Poisson's ratio for case 1 are:  1 = 250 MPa,  2 = 2.5 MPa,  1 = 0.5 and  2 = 0.5.For case 2,  2 = 0.2 with other parameters being the same as those of case 1.The load value is  = 1 N/mm for both cases.
The convergence of the Y-displacement of point A with respect to mesh refinement, as shown in Fig. 19, follows the same trend for both cases.Interestingly, the P2/P1 and P3/P1dc elements converge from above, while Q1/P0 and P2b/P1dc elements converge from below.Although the displacement solution obtained with all elements converges with mesh refinement, the pressure field obtained with the P2/P1 element differs significantly at the interface due to the continuity of the pressure DOFs, see Fig. 20.Even though P2b/P1dc and P3/P1dc elements capture the pressure jump accurately, the pressure field obtained with the P3/P1dc element is smoother within each domain.

Block with a single material
We first consider the benchmark example of a homogeneous block under compression.Due to the symmetry of the boundary and loading conditions, a quarter portion of the block with symmetry boundary conditions is modelled, as shown in Fig. 21(a).The material model is the Neo-Hookean model with Young's modulus of 240.565MPa and a Poisson's ratio of 0.4999.The problem is studied using five meshes shown in Fig. 22, and the Z-displacement of point A as the percentage of compression is plotted in Fig. 21(b) for different elements.As shown, accurate results are obtained with coarse meshes using P2b/P1dc and P3/P1dc elements.From the contour plots of pressure presented in Fig. 23, we can see that the pressure field obtained with the P3/P1dc element is smoother than that of the P2b/P1dc element.
Computer Methods in Applied Mechanics and Engineering 432 (2024) 117354

Block with four layers
Next, we consider the block with four layers from Wu et al. [12].The setup of the problem is shown in Fig. 24.The material properties of each layer are tabulated in Table 1.The analysis is performed using two meshes shown in Fig. 24 with P2/P1, P2b/P1dc and P3/P1dc elements.The first and second meshes consist of 627 and 3742 elements, respectively.The number of nodes for the respective meshes are 1131 and 5935 for the P2 element, 1758 and 9677 for the P2b element, and 3475 and 18935 for the P3 element.The load is applied in five uniform increments.The deformed shapes are presented in Fig. 25.The graphs of Z-displacement along the diagonal like AC presented in Fig. 26 show that the results obtained with different elements match well with the reference solution obtained with finer meshes (48004 nodes) in Wu et al. [12].The contour plots of the pressure field in Fig. 27 show that P2b/P1dc and P3/P1dc elements capture the pressure jumps accurately while the P2/P1 element fails to do so.

Twisting of a layered column
We now consider the torsion of a square column with four layers.The problem setup and the meshes used for the analysis are presented in Fig. 28.The material model is the incompressible Neo-Hookean with  = 170 MPa for layers 1 and 3 and  = 17 MPa for layers 2 and 4. The simulations are carried out for one full rotation with a load step of 9 degrees.The deformed shapes and contour plots of pressure obtained with the different element types at three different load steps are presented in Fig. 29.While there is a negligible difference in the deformed shapes, the pressure field obtained with the P3/P1dc element is smoother within the individual domains when compared with the P2b/P1dc element.
The wall times to complete 40 load steps with P2/P1, P2b/P1dc and P3/P1dc elements are, respectively, 75, 87 and 510 s using the MUMPS direct solver on a single processor.It is understandable that P3/P1dc elements are computationally more expensive than the P2/P1 and P2b/P1dc elements; however, the accuracy of results is substantially higher.P2b/P1dc element does offer a compromise between the accuracy and computational cost while accurately maintaining the pressure and stress jumps across the   material interfaces; however, the accuracy of pressure fields obtained within individual domains is subpar when compared with the P3/P1dc element.

Twisting of a fibre-reinforced soft tissue
We now consider a soft tissue sample of length 5 mm and diameter 2 mm that is reinforced with fibres of diameter 0.4 mm placed at (±0.3, ±0.3) mm, as depicted in Fig. 30.The material model is the incompressible Neo-Hookean model with the shear modulus of 1 MPa and 10 MPa, respectively, for the matrix and fibre.The problem is simulated with Q1/P0, P2b/P1dc and P3/P1dc elements using the meshes shown in Fig. 30 for a rotation of 360 degrees using increments of 9 degrees.The solution obtained with the Q1/P0 element is used as the reference.The deformed shapes at 180 and 360 degrees obtained with the Q1/P0 and P3/P1dc elements are shown in Fig. 31.The contour plots of pressure presented in Figs.32 and 33 clearly show that the pressure field obtained with the P2b/P1dc element is significantly poor when compared with the P3/P1dc element.

Particle-filled soft composite
As the last example, we consider a soft matrix filled with particles.This example demonstrates the need for accurate and robust tetrahedral elements since the generation of meshes with hexahedral elements for such problems is cumbersome.The matrix is a cube of side length 1 mm and is filled with six spherical particles of radii 0.1, 0.12, 0.1, 0.12, 0.1 and 0.12 mm with their respective centres at (0.25,0.6,0.1),(0.6, 0.5, 0.3), (0.6, 0.3, 0.7), (0.3, 0.3, 0.3), (0.3, 0.7, 0.3) and (0.7, 0.7, 0.7).The material model for the matrix is incompressible Mooney-Rivlin with  1 = 0.375 and  2 = 0.125 MPa, and that for the particles is incompressible Neo-Hookean with  = 100 MPa.The face at  = 0 is fixed, and a displacement of (0.8, 0, 0) mm is applied on the face at  = 1 mm in increments of 0.1 mm.The mesh used for the simulation is shown in Fig. 34.P3/P1dc element are much smoother than those of the P2b/P1dc element, indicating better accuracy observed in the thick-walled cylinder and sphere examples with linear elasticity.

Conclusions
We have assessed the performance of discontinuous pressure spaces based on triangular and tetrahedral elements in the context of mixed displacement-pressure formulation for the simulation of multimaterial problems in linear and nonlinear elasticity.The novel contribution of the present work is the adaptation of P2b/P1dc and P3/P1dc elements and extensive studies of their relative advantages and disadvantages for problems in linear and nonlinear elasticity.The accuracy of different elements in accurately resolving pressure jumps across material interfaces, as well as in capturing pressure fields within the individual material domains, is assessed using several numerical examples.
The key observations from the results obtained in the present work are: • The P2/P1 element, without any additional preprocessing, fails to capture pressure jumps across material interfaces due to the continuity of pressure nodes/DOFs across element edges/faces.• The magnitudes of error norms in pressure obtained with the P2b/P1dc element are higher than those obtained with the Q1/P0 element despite the P2b/P1dc element having significantly higher pressure DOFs.• For the single material case, converting the P2/P1 element to the P2b/P1dc element by adding an extra node with a cubic bubble function does not lower the displacement error norm significantly.However, surprisingly, the P2b/P1dc element yields pressure fields of substantially poor quality despite a significant increase in pressure DOFs.• Even though the P2b/P1dc element resolves the pressure jumps across material interfaces accurately, the accuracy of pressure fields within the individual domain is poor.Despite having the same number of pressure DOFs, the magnitudes of error in  pressure obtained with the P2b/P1dc element are about two orders of magnitude higher when compared with the P3/P1dc element.The superior accuracy of the P3/P1dc element is evident in the smoother pressure fields obtained with it, as shown with the warped plots of the pressure field for the thick-walled cylinder example.These differences in the accuracy have a substantial effect on the pressure fields obtained for the nonlinear elasticity problems, as demonstrated by the examples of twisting of a layered column and fibre-reinforced soft tissue.• While the computational cost of the P3/P1dc element is higher than that of the P2b/P1dc element, the accuracy of results obtained with P3/P1dc is substantially better.
In conclusion, the present work offers significant insights into the performance of Q1/P0, P2/P1, P2b/P1dc and P3/P1dc elements for simulating multimaterial problems in linear and nonlinear elasticity.While the Q1/P0 element is reasonably accurate and competitive in terms of performance, the difficulty of mesh hexahedral generation for complex problems makes it less appealing for practical problems.The P2b/P1dc element performs poorly, requiring finer meshes for accurate results.The P3/P1dc yields accurate results using coarse meshes.We hope this work inspires further research on new element technologies for multimaterial problems consisting of incompressible hyperelastic materials.
3(b), to demonstrate the effectiveness of different element types in capturing the pressure fields and stress discontinuities at the material interfaces.The radii of the circular portions are   = 100 mm,  0 = 200 mm and   = 150 mm.The material properties are:  = 200 MPa,  = {0.3,0.5}, and  1 = 200 MPa,  2 = 20 MPa and  1 =  2 = {0.3,0.5}.The pressure on the inner surface of the cylinder is   = 0.1 MPa.

Fig. 14 .
Fig. 14.Thick-walled composite sphere: contours of plots of pressure obtained with different element types using mesh 2.

Fig. 15 .
Fig. 15.Cook's membrane: setup of (a) single material case and (b) two material case.

Fig. 16 .
Fig. 16.Cook's membrane with a single material: meshes used for the analysis.

ComputerFig. 17 .
Fig. 17.Cook's membrane with a single material: convergence of Y-displacement of point A for  = 240.565MPa and  = 0.4999 with (a) linear elastic material and Neo-Hookean material.Reference values are from Kadapa et al. [15].

Fig. 18 .
Fig. 18.Cook's membrane with a single material: contours of pressure obtained with the Neo-Hookean model and different element types for  = 0.4999.

Fig. 20 .
Fig. 20.Cook's membrane with two materials: contours of pressure obtained with the linear elastic material model and different element types.

Fig. 23 .
Fig. 23.Block under compression: contour plots of pressure obtained with different element types for the mesh with 12 elements per side.

Fig. 25 .
Fig. 25.Multilayer block: deformed shapes and contour plots of the Z component of displacement obtained with different elements.

Fig. 27 .
Fig. 27.Multilayer block: contour plots of pressure obtained with different element types with mesh 2.

Fig. 28 .
Fig. 28.Twisting of a layered column: (a) problem setup and (b) the mesh used for the simulation.

Fig. 33 .
Fig. 33.Twisting of a fibre-reinforced soft tissue: contour plots of pressure for the clipped domains.

Fig. 36 .
Fig. 36.Particle-filled composite: contour plots of pressure for the particles at 0.2 mm and 0.8 mm of applied displacement.

Fig. 37 .
Fig. 37. Particle-filled composite: contour plots of pressure on the clipped model (with a plane at 0.4 mm along the Y-axis) at 0.2 mm and 0.8 mm of applied displacement.