Within this paper we discuss a numerical strategy to solve the elasticity problem upon unstructured and non conforming meshes, allowing all kinds of flat-faced elements (polygons in 2D and polyhedra in 3D). The core of the formulation relies on two numerical procedures 1) the Control Volume Function Approximation (CVFA), and 2) the polynomial interpolation in the neighborhood of the control volumes, which is used to solve the surface integrals resulting from applying the divergence theorem. By comparing the estimated stress against the analytical stress field of the well known test of an infinite plate with a hole, we show that this conservative approach is robust and accurate.
Keywords: Stress analysis, elasticity, solid mechanics, CVFA, Finite Volume, unstructured mesh, non conforming mesh
One of the main aims in engineering is creating tools, structures and systems to enhance the quality of life in our society. In the course of the creation process, the design stage is critical for the final outcome. During this stage the engineer have to predict the prototype response when interacting with the physical world. Many of the observed phenomena in the physical world, such as solid mechanics, fluid dynamics, heat diffusion, and others, can be described with Partial Differential Equations (PDEs) by assuming time and space as a continuum.
Computational Continuum Mechanics (CCM) is the area dedicated to develop numerical methods and heuristics to solve these PDEs. Most of the methods can be classified into these two families: 1) weighted residual and 2) conservative methods. The Galerkin formulations are popular and widely used weighted residual methods, such as the Finite Element Method (FEM), which is a well established technique in Computational Solid Mechanics (CSM). Alternatively, the Finite Volume Method (FV) and the Control Volume Function Approximation (CVFA) are common approaches of conservative methods. The main difference between both families is that weighted residuals methods do not conserve quantities locally, but globally instead, resulting in linear systems with commendable numerical properties (symmetrical and well-conditioned matrices, for example). Nevertheless, due to its conservative nature, the second group is more attractive for fluid structure interaction [1,2] and multiphysics simulations [3,4], where several PDE-solvers must be coupled. For that reason, in recent years FV has been subject of interest for solving CSM problems.
Most of the CSM non-linear strategies depend on the accuracy of the estimated stress field for the elasticity problem, such as those for plasticity and damage [5,6]. Hereafter we refer as elasticity-solver to the numerical computation that calculates the displacement and stress fields for a given domain and boundary conditions.
In his influential papers, Oñate et al. propose a FV format for structural mechanics based on triangular meshes [7,8], discussing the cell vertex scheme, the cell centred finite volume scheme and its corresponding mixed formulations, showing that the cell centred strategy produces the same symmetrical global stiffness matrix that FEM using linear triangular elements. Analogously, Bailey et al. develop a similar approach, but using quadrilateral elements to produce cell-centred volumes [9,10]. Even though, the shapes of the volumes in both formulations are completely defined by the FEM-like mesh (triangular or quadrilateral) and it is not possible to handle arbitrary polygonal shapes, as we might expect when receiving the mesh as a read-only input.
Slone et al.  extends the investigation of  by developing a dynamic solver based on an implicit Newmark scheme for the temporal discretization, with the motivation of coupling an elasticity-solver with his multi-physics modelling software framework, for later application to fluid structure interaction.
Another remarkable algorithm is the proposed by Demirdzic et al. [12,13,14,15,16]. The numerical procedure consists in decoupling the strain term into the displacement Jacobian and its transpose in a cell-centred scheme. The Jacobian is implicitly estimated by approximating the normal component of each face as the finite difference with respect to the adjacent nodes, while the Jacobian transpose is an explicit average of Taylor approximations around the same adjacent nodes. This decoupling produces a smaller memory footprint than FEM because the stiffness matrix is the same for all the components. The solution is found by solving one component each iteration in a coordinate descent minimization. This line of work has shaped most of the state of the art techniques in FV for coupling elasticity-solvers to Computational Fluid Dynamics (CFD) via finite volume practices (usually associated to CFD), such as the schemes proposed by [17,18,19]. However, this segregated algorithm may lead to slow convergence rates when processing non-linear formulations, for example, when it is required to remove the positive principal components of the stress tensor in phase-field damage formulations . In addition, if some non-linear strategy requires multiple iterations of the linear elasticity-solver, such as finite increments in damage models, the nested iterations will increase the processing requirements for simple problems. To circumvent this drawback Cardiff et al. presents a fully block coupled direct solution procedure , which does not require multiple iterations, at expense of decomposing the displacement Jacobian of any arbitrary face into a) the Jacobian of the displacement normal component, b) the Jacobian of the displacement tangential component, c) the tangential derivative of the displacement normal component and d) the tangential derivative of the displacement tangential component. This decomposition complicates the treatment of the stress tensor in the iterative non-linear solvers mentioned before for plasticity and damage.
A generalized finite volume framework for elasticity problems on rectangular domains is proposed by Cavalcante et al. . They use higher order displacement approximations at the expense of fixed axis-aligned grids for discretization.
Nordbotten proposes a generalization of the multi-point flux approximation (MPFA), which he names multi-point stress approximation (MPSA) . The MPSA assembles unique linear expressions for the face average stress with more than two points in order to capture the tangential derivatives. The stress field calculated with this procedure is piece-wise constant.
In this work, we propose an elasticity-solver based on CVFA techniques [24,25], using piece-wise polynomial interpolators for solving the surface integrals on the volumes boundaries. The polynomials degree can be increased without incrementing the system degrees of freedom, which make this method more suitable for non-linear models and dynamic computations. Furthermore, this algorithm can handle polygonal/polyhedral, unstructured and non conforming meshes, and does not require the decomposition of the stress tensor.
We consider an arbitrary body, , with boundary . The displacement of a point is denoted by . The subscript brackets indicates the component of the vector or matrix. By assuming small deformations and deformation gradients, the infinitesimal strain tensor, denoted , is given by
Moreover, by considering isotropic elastic materials, the stress tensor, , is defined as
where is the identity matrix, is the trace, and and are the Lamé parameters characterizing the material. These parameters are related with the Young's modulus, , and the Poisson's ratio, , by the following equivalences
where for plane stress analysis, and for plane strain and the 3D cases.
The strong form of the static linear elasticity equation is
The domain boundary is the union of the boundary with Neumann conditions, denoted , and the boundary with Dirichlet conditions, denoted . The equation (5) must be satisfied for the given forces and displacements along the boundary. The following equations remark these user defined boundary conditions,
The Figure 1 illustrates the initial body with the boundary conditions, and the distorted body after equilibrium is solved in the elasticity equation.
|Figure 1. (a) The initial body, , with its boundary conditions . (b) The distorted body resulting from solving equilibrium in the elasticity equation, with the boundary conditions given by and . The boundaries where there are not conditions indicated explicitly, correspond to Neumann conditions|
On this section we go into the details of the numerical procedure by discussing 1) the discretization with CVFA, 2) the control volumes integration, 3) the subfaces integrals, 4) the simplex-wise polynomial approximation, 5) the pair-wise polynomial approximation, 6) the homeostatic splines used within the shape functions, 7) the linear system assembling, 8) how to impose boundary conditions and 9) two special cases of the formulation.
For the sake of legibility, in some parts of the text we unfold the matrices only for the bidimensional case, but the very same procedures must be followed for the 3D case.
The domain is discretized into control volumes, denoted , using the Control Volume Function Approximation (CVFA) proposed by Li et al. [24,25]. The partition of is defined by
where the boundary of each control volume, , is composed by flat faces, denoted ,
|Figure 2. The partition is the discretization of the domain into control volumes. The boundary of the control volumes, , is conformed by flat faces, denoted . The unit vector is normal to the face . The faces of the volumes adjacent to the boundary are integrated using the condition|
The Figure 3 shows a three dimensional control volume.
|Figure 3. The boundary of the three dimensional control volume is subdivided into flat faces, denoted . The unit vector is normal to the face|
Every control volume must have a calculation point
which is used to estimate the displacement field. Such a point is the base location to calculate the stiffness of the volume. In the volumes adjacent to the boundary , it is convenient to establish the calculation point over the corresponding boundary face,
in order to set the Dirichlet condition directly on the point.
The elasticity equation (5) is now integrated over the control volume
using the divergence theorem we transform the volume integral into a surface integral over the volume boundary
The evaluation of the surface integrals is based on the approximation of the displacement field inside the neighborhood of the volume, denoted ,
making use of a group of piece-wise polynomial interpolators, denoted . We are going to discuss these interpolators later in this section.
For that reason, the displacement field is decoupled from the stress tensor by using the strain (1) and stress (2) definitions. Taking advantage of the stress tensor symmetry , we rewrite the stress normal to the boundary as
where is the face orientation matrix and is the engineering stress vector. Developing the stress definition (2) component-wise we can decompose it into the constitutive matrix, denoted , and the engineering strain vector, denoted , as follows
then the components of the strain vector are retrieved from the equation (1), and it is decomposed into the matrix differential operator and the displacement function .
where is the stiffness of the volume boundary.
Once the displacement field is decoupled, we rewrite the equation (12) as
Notice that the face orientation along the flat face, denoted , is constant. Furthermore, if the control volumes are considered to be made of a unique material and the flat faces to be formed by pairs of adjacent volumes, then the constitutive matrix along the flat face, denoted , is also considered constant. The matrix is estimated from the harmonic average of the Lamé parameters assigned to the adjacent volumes, where and are the properties of the volume ,
with and we simplify the equation (20) as
where is the strain integral along the flat face ,
The accuracy of the method depends on the correct evaluation of this integral.
The surface integrals along the flat faces are calculated using an auxiliary piece-wise polynomial approximation of the displacement field. This approximation is based on the simplices (triangles in 2D or tetrahedra in 3D) resulting from the Delaunay triangulation of the calculation points from the neighborhood of . The Delaunay triangulation is the best triangulation for numerical interpolation, since it maximizes the minimum angle of the simplices, which means that its quality is maximized as well. We define the neighborhood of volume as the minimum set of calculation points such that the simplices intersecting does not change if we add another calculation point to the set. Observe that the neighborhood does not always coincide with the set of calculation points in volumes adjacent to , as in most of the FV formulations. Once the neighborhood is triangulated, we ignore the simplices with angles smaller than degrees, and the simplices formed outside the domain, which commonly appear in concavities of the boundary . The local set of simplices resulting from the neighborhood of is denoted . The Figure 4 illustrates the difference between (a) the simplices resulting from the triangulation of the calculation points in adjacent volumes and (b) those resulting from the triangulation of the proposed neighborhood .
|Figure 4. (a) The dotted line illustrates the triangulation of the calculation points of adjacent volumes to , used by most of the FV methods. (b) The dotted line shows the simplices forming the piece-wise approximation used to solve the integrals of the control volume|
The face is subdivided into subfaces, denoted ,
these subfaces result from the intersection between and the control volume . The Figure 4.b illustrates six key points of this approach, 1) the simplices are used to create a polynomial interpolation of over the boundary of the control volume, 2) most of the faces are intersected by several simplices, such faces must be divided into subfaces to be integrated, 3) some few faces are inside a single simplex, as illustrated in the face formed by and , 4) there are volumes that require information of non-adjacent volumes to calculate its face integrals, such as requires , 5) the dependency between volumes is not always symmetric, which means that if requires does not implies that requires , and 6) non conforming meshes are supported, as shown in the faces formed by , , , and .
The integral (23) is now rewritten in terms of the subfaces
Each subface is bounded by a simplex, where the displacement , and it derivatives, , are a polynomial interpolation. Hence the integrals in equation (25) are solved exactly by using the Gauss-Legendre quadrature with the required number of integration points, denoted , depending on the polynomial degree,
|Figure 5. (a) The shaded volume is being integrated. The integral over the subface is calculated using the polynomial approximation of the shaded simplex. The integration point must be mapped to (b) the normalized space in order to use the Gauss- Legendre quadrature|
Most of the cases, the displacement is interpolated inside the simplices, but in some geometrical locations these can not be created, in consequence, the displacement is interpolated pair-wise using the volumes adjacent to the subface . We discuss both strategies in the following subsections.
In the general case, the simplices are formed by points. The points forming the simplex that is bounding the subface are denoted , and its displacements .
The shape functions used for the polynomial interpolation are defined into the normalized space. A point in such space is denoted , its component is denoted , and the point forming the simplex is expressed . The nodes of the normalized simplex are given by the origin, , and the standard basis vectors,
|Figure 6. (a) The simplex formed by the points , and in the original space contains an interior point that is mapped to (b) into the normalized 2D-simplex formed by the points , and|
|Figure 7. (a) The 3D-simplex formed by the points , , and in the original space contains an interior point that is mapped to (b) into the normalized 3D-simplex formed by the points , , and|
The shape functions, denoted , are used to interpolate the displacement field inside the normalized simplex. Such functions are non-negative and are given by
where is the homeostatic spline, which is the simplest polynomial defined in the interval that have derivatives equal to zero in the endpoints of the interval. We will discuss this spline later.
The set of shape functions is a partition of unity, which means that the sum of the functions in the set is equal to one into the interpolated domain
furthermore, these functions are equal to one in its corresponding node, which implies that
The gradients of the shape functions with respect to the normalized space are denoted . The norm of the sum of such gradients is zero
which means that there are not numerical artifacts into the strain field.
Any point inside the simplex can be formulated as a function of a point in the normalized space, , by using the shape functions and the points forming the simplex
In order to calculate the normalized point, denoted , associated to the integration point , we use the shape functions definitions to rewrite the equation (33) in matrix form
where is the vector resulting from evaluating the spline for component-wise, and is the distortion matrix given by the concatenation of the following column vector differences
Now, from equation (35) we retrieve the point as
and solving for we obtain
where is the inverse function of applied component-wise to the product of the matrix-vector operation.
Similar to the approximation in equation (33), within the simplex enclosing the subface , the displacement field evaluated at is defined as,
Hence, when calculating the quadrature of equation (26), the strain evaluated at the integration point is given by
where captures the deformation at , and is the vector with the concatenated displacement components of the points forming the simplex.
In order to calculate the deformation matrix , we require the derivatives of the shape functions with respect to , denoted . These derivatives are calculated by solving the linear systems resulting from the chain rule
where is the geometric jacobian evaluated at . This jacobian relates both spaces, captures the distortion of the simplex, and is derivated from equation (33),
The gradients of the shape functions with respect to inside the sum are obtained straightforward once we have the spline first derivative . Notice that the geometric jacobian is equivalent to the distortion matrix if and only if the homeostatic spline is .
Since we are not making any assumption about the volumes distribution through the mesh, neither about the internal location of its calculation points, then we have to deal with portions of the mesh that are no covered by any simplex. The Figure 8 illustrates the two most common cases. The first case takes place in meshes where the calculation points of volumes contiguous to the boundary are in the interior of such volumes, producing subfaces not intersected by any simplex, and the second case occurs when elongated sections of the domain are discretized with a queue of aligned volumes, where each volume has only two neighbors on opposite faces and no simplex can be formed.
|Figure 8. (a) When the calculation points of volumes contiguous to the boundary are in the interior of such volumes, there will arise subfaces next to the boundary that can not be covered by any simplex. (b) Portions of the mesh formed by a queue of aligned volumes do not allow the formation of simplices through that queue and there will be whole faces not covered by any simplex|
In such cases, the displacement field within the subface is a pair-wise polynomial approximation between the adjacent volumes, and , regardless the dimension
where and are the shape functions, and is the normalized projection of the integration point into the vector which goes from to , denoted
When calculating the quadrature of equation (26), the pairwise strain is given by
|Figure 9. The gradient of the pairwise approximation is not constant along the face , since its normal is not necessary aligned with . The integration point is projected into to evaluate the deformation matrix|
This pairwise approximation must be used only when necessary because it can not capture the deformation orthogonal to .
The homeostatic spline is a function of a single variable defined from to , denoted , and curved by the parameter , which indicates the level of smoothness. This spline is the simplest polynomial with derivatives equal to zero at the endpoints and . The polynomial degree is given by , and such a polynomial requires integration points to calculate the exact integral in equation (26) using the Gauss-Legendre quadrature.
When designing this spline, we wanted to gain accuracy by building a piece-wise bell-shaped interpolation function around the calculation points, inspired on the infinitely smooth kernels used in other numerical techniques. Therefore, we force the derivatives of the polynomial to be zero over such points in order to homogenize the function. For that reason, we use the term homeostatic spline when referring to this spline.
To fulfill the smoothness requisites commented before, we solved a linear system for calculating the coefficients of the polynomial. The equations of this system were obtained by forcing the derivatives to be zero at the endpoints. Once we solved the coefficients for the first twenty polynomials, from to , we found out that the first half of such coefficients are null, and the entire polynomial can be calculated directly as
where is the not null coefficient
is the number of factors needed to calculate
and is accumulation of the coefficients for normalizing the spline
The first derivative is simply calculated as
|Figure 10. (a) The evolution of the homeostatic spline from to illustrates the smoothness requirements at the endpoints of each spline and its (b) first derivatives|
Since the derivatives of the homeostatic spline 50 are zero at the endpoints of the interval , the inverse function is not defined in that points. However, we estimate a pseudo-inverse within this interval, , by finding the coefficients of a polynomial of the same degree, , such that the endpoints coincide with the spline and the first derivative at the midpoint is equivalent to the inverse of the spline first derivative, that is
The higher derivatives in the midpoint are forced to be zero. Once we calculated the coefficients for the first twenty polynomials, from to , we found out that the pseudo-inverse can be approximated directly from the following formulae
where is the coefficient for , and is the factor that distinguish higher order coefficients. Such terms are calculated as
respectively. The Figure 11 exhibits the curves for the first seven levels of smoothness. The null higher derivatives requirement is noticeable at the midpoint.
|Figure 11. Curves of the pseudo-inverse for the first seven levels of smoothness. The slope at the midpoint exposes the null higher derivatives requirement when increasing the polynomial order|
The Figure 12 shows the shape functions for the 2D case. The top displays the last node function and the bottom the first node function, the function of the second node is equivalent to that of the first one. The columns separate the first three levels of smoothness. Top and bottom functions coincides at the edges in order to create a continuous field, but only the bottom functions decay uniformly from the node to the opposite edge. The shape functions with are the only case where all the shape functions are indistinguishable, these are planes.
|Figure 12. For the bidimensional case, the top displays the last node function and the bottom the first node function, the function of the second node is equivalent to that of the first one. The columns separate the first three levels of smoothness|
The Figure 13 shows the magnitude of the gradient with respect to the normalized space. With the same tabular configuration of Figure 12, the columns separate the first three levels of smoothness, the top displays the last node gradient and the bottom the first node gradient, the gradient of the second node is equivalent to that of the first one. Only the gradient magnitude at the bottom has a uniform variation from the node to the opposite face, and the value of the node does not contribute to the value of such a face. On the contrary, in the top can be observed that the value of the node contributes to the gradient at the opposite face, which means that using the continuity on the stress field is only guaranteed at the calculation points, but not in the simplices edges.
|Figure 13. For the bidimensional case, the top displays the last node gradient magnitudes and the bottom the first node gradient magnitudes, the gradient magnitudes of the second node is equivalent to that of the first one. The columns separate the first three levels of smoothness|
then, the volume equilibrium equation (22) is
reordering the terms we obtain
where the matrix
is the stiffness contribution at the integration point , within the subface when integrating the volume. Observe that the stiffness matrix is rectangular and the resulting global stiffness matrix is asymmetric.
The Neumann boundary conditions are imposed over the volume faces intersecting the boundary, by replacing the corresponding term in the sum of equation (20) with the integral of the function provided in (6.a),
The Dirichlet conditions are imposed over the volumes calculation points by fixing the displacement as it is evaluated in the function given in (6.b),
Since the Dirichlet conditions are imposed directly on the calculation points, these points must be located along the face which intersects the boundary with the condition .
By making some considerations, we identify two special cases where the calculations can be simplified, in order to increase the performance of the total computation, at the expense of losing control over the volumes shape. These cases are 1) the Voronoi mesh assumption and 2) the FV-FEM correlation.
In the first case, we assume that the initial mesh is equivalent to the Voronoi diagram and that the Voronoi centres correspond to the calculation points . Hence, the subdivision of the neighborhood is already given by the Delaunay triangulation which is dual to the Voronoi mesh, as illustrated in the Figure 14.a. Moreover, the integrals of subfaces using pair-wise approximations can be exactly integrated with the midpoint rule, since the faces are orthogonal to the vector joining the calculation points , and the derivatives along the subface are constants.
|Figure 14. (a) The initial mesh is equivalent to the Voronoi diagram and the Voronoi centres correspond to the calculation points . (b) The initial mesh is generated from a FEM-like triangular mesh. The calculation points are defined to be the nodes of the triangular mesh, and the volume faces are created by joining the centroids of the triangles with the midpoint of the segments|
In the second case, the initial mesh is generated from a FEM-like triangular mesh and the approximations are assumed to be linear. In such a case, the calculation points are defined to be the nodes of the triangular mesh, and the volume faces are created by joining the centroids of the triangles with the midpoint of the segments, as presented in Figure 14.b. This particular version is equivalent to the cell-centred finite volume scheme introduced by Oñate et al. , who proved that the global linear system produced by this FV scheme is identical to that produced by FEM if the same mesh is used.
In order to test the numerical performance of the proposed method, we use the well known analytical experiment of an infinite plate with a circular hole in the origin . In such a experiment, the plate is stretched along the horizontal axis with a uniform tension of from each side, as is shown in Figure 15. The material is characterized by the Poisson ratio, , and Young modulus, . Plane stress is assumed with thickness equivalent to the unity. The dimensions of the computational domain are and .
|Figure 15. (a) Infinite plate with a hole being stretched along the horizontal axis with a force of from each side. (b) Computational domain, and , with axysymmetrical assumptions used to test the numerical method. The polar coordinates, and , for calculating the analytical stress field|
The analytical solution is given by the following formulae
where the polar coordinates,
are used within the calculus. The Figure 16 exhibits (a) the discretization of the computational domain into 2411 polygonal volumes used to compare the numerical results against the analytical stress field. This mesh is not equivalent to the Voronoi diagram. (b) Level sets of between to , with steps of . (c) Level sets of between and , with steps of . (d) Level sets of between and , with steps of .
|Figure 16. (a) Polygonal mesh used for comparison of numerical results. (b) Level sets of between to . (c) Level sets of between and . (d) Level sets of between and|
The Dirichlet conditions are imposed on the bottom and left side of the computational domain as is shown in the Figure 16.b. Next in order, the analytic stress of equations (64), (65) and (66) is imposed as Neumann condition over the top and right side of the computational domain.
The Figure 17.a presents the averaged error as a function of the volumes face length mean, denoted , as we might expect, the error is proportional to the mesh refinement. For a mesh of 478 volumes, the Figure 17.b shows the percentage of the error with respect to the error of , for different smoothing levels, correspond the linear interpolator. Observe that the error of the stress field does not decreases significantly because we do not increase the degrees of freedom of the linear system, although we built a better field description, which can be useful when solving non-linear formulations.
|Figure 17. (a) The averaged error decreasing as a function of mesh size, denoted . (b) Using a mesh of 628 volumes, percentage of error for different smoothing levels with respect to the error of , which is the linear interpolator, error increases after|
In this work we discussed a numerical technique to solve the elasticity equation for unstructured and non conforming meshes formed by elements of any arbitrary polygonal/polyhedral shape. The elasticity-solver is based on a finite volume formulation that, using the divergence theorem, represent the volume integral of the stress divergence in terms of the surface integral of the stress over the volume boundary. For considering the unit vector normal to the boundary as a constant, the boundary is divided into flat faces. Conforming and non-conforming meshes are processed without distinction. The displacement field is a piece-wise polynomial approximation surrounding the volumes, built on the top of the simplices resulting from the Delaunay triangulation of the volume neighborhood. A pair-wise polynomial interpolation is used for neighborhoods where the simplices are not the best option, and for 1D problems.
We introduced the homeostatic splines and it pseudo-inverse for higher order polynomial interpolations without the necessity of increasing the discretization points.
Since the stress term is calculated directly on the boundary of the control volumes, this strategy can be used in fracture formulations where the volumes are treated as indivisible components and the rupture occurs across the volumes boundaries.
In future work, we will investigate the accuracy of this numerical procedure for solving non-linear and dynamic models when using higher order homeostatic splines.
The present investigation was sponsored by a CONACYT scolarship from the Mexican government and the TCAiNMaND project, an IRSES Marie Curie initiative under the European Union 7th Framework Programme. In addition, the authors want to express his gratitude to MSc. Ernesto Ortega, and to Drs. Rafel Herrera and Joaquin Peña for his insightful commments and discussions during the elaboration of this paper.
 Slone A.K., Pericleous K., Bailey C., Cross M. Dynamic fluid-structure interaction using finite volume unstructured mesh procedures. Computers and Structures, 80:371-390, 2002.
 Slone A.K., Pericleous K., Bailey C., Cross M., Bennett C. A finite volume unstructured mesh approach to dynamic fluid-structure interaction: an assessment of the challenge of predicting the onset flutter. Applied Mathematical Modeling, 28:211-239, 2004.
 Bailey C., Taylor G.A., Cross M., Chow P. Discretisation procedures for multi-physics phenomena. Journal of Computational and Applied Mathematics, 103:3-17, 1999.
 Souli M., Ouahsine A., Lewin L. ALE formulation for fluid-structure interaction problems. Computer Methods in Applied Mechanics and Engineering, 190:659-675, 2000.
 Lubliner J., Oliver J., Oller S., Oñate E. A Plastic-Damage model for concrete. International Journal of Solids and Structures, 25:299-326, 1989.
 Armero F., Oller S. A general framework for continuum damage models. I. Infinitesimal plastic damage models in stress space. International Journal of Solids and Structures, 37:7409-7436, 2000.
 Oñate E., Cervera M., Zienkiewicz O.C. A finite volume format for structural mechanics. International Journal for Numerical Methods in Engineering, 37:181-201, 1994.
 Idelsohn S.R., Oñate E. Finite volumes and finite elements: Two 'Good Friends'. International Journal for Numerical Methods in Engineering, 37:3323-3341, 1994.
 Fryer Y.D., Bailey C., Cross M., Lai C.H. A control volume procedure for solving the elastic stress-strain equations on an unstructured mesh. Applied Mathematical Modeling, 15:639-645, 1991.
 Bailey C., Cross M. A Finite Volume procedure to solve elastic solid mechanics problems in three dimensions on an unstructured mesh. International Journal for Numerical Methods in Engineering, 38:1757-1776, 1995.
 Slone A.K., Bailey C., Cross M. Dynamic solid mechanics using finite volume methods. Applied Mathematical Modeling, 27(2):69-87, 2003.
 Demirdzic I., Martinovic D. Finite volume method for thermo-elasto-plastic stress analysis. Computer Methods in Applied Mechanics and Engineering, 109:331-349, 1993.
 Demirdzic I., Muzaferija S. Finite volume methods for stress analysis in complex domains. International Journal for Numerical Methods in Engineering, 137:3751-3766, 1994.
 Demirdzic I., Muzaferija S., Peric M. Benchmark solutions of some structural analysis problems using finite-volume method and multigrid acceleration. International Journal for Numerical Methods in Engineering, 40:1893-1908, 1997.
 Bijelonja I., Demirdzic I., Muzaferija S. A finite volume method for large strain analysis of incompressible hyperelastic materials. International Journal for Numerical Methods in Engineering, 64:1594-1609, 2005.
 Bijelonja I., Demirdzic I., Muzaferija S. A finite volume method for incompressible linear elasticity. Computer Methods in Applied Mechanics and Engineering, 195:6378-6390, 2006.
 Jasak H., Weller H.G. Application of the finite volume method and unstructured meshes to linear elasticity. International Journal for Numerical Methods in Engineering, 48:267-287, 2000.
 Tukovic Z., Ivankovic A., Karac A. Finite-volume stress analysis in multi-material linear elastic body. International Journal for Numerical Methods in Engineering, 93(4):400-419, 2012.
 Tang T., Hededal O., Cardiff P. On finite volume method implementation of poro-elasto-plasticity soil model. International Journal for Numerical and Analytical Methods in Geomecanics, 39:1410-1430, 2015.
 Borden M.J., Verhoosel C.V., Scott M., Hughes T.J.R., Landis C.M. A phase-field description of dynamic brittle fracture. Computer Methods in Applied Mechanics and Engineering, 217-220:77-95, 2012.
 Cardiff P., Tukovic Z., Jasak H., Ivankovic A. A block-coupled Finite Volume methodology for linear elasticity and unstructured meshes. Computers and Structures, 175:100-122, 2016.
 Cavalcante M.A.A., Pindera M-J. Generalized Finite-Volume Theory for Elastic Stress Analysis in Solid Mechanics - Part I: Framework. Journal of Applied Mechanics, 79(5):1006-1017, 2012.
 Nordbotten J.M. Cell-centered finite volume discretizations for deformable porous media. International Journal for Numerical Methods in Engineering, 100:399-418, 2014.
 Li B., Chen Z., Huan G. Control volume function approximation methods and their applications to modeling porous media flow. Advances in Water Resources, 26(4):435-444, 2003.
 Li B., Chen Z., Huan G. Control volume function approximation methods and their applications to modeling porous media flow II: the black oil model. Advances in Water Resources, 27(2):99-120, 2004.
 Timoshenko S.P., Goodier J.N. Theory of Elasticity. McGraw-Hill, 3 Edition, 1970.