Finite volume for stress analysis upon polygonal, unstructured and non conforming meshes

==Abstract==

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

==1. Introduction==

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 <span id='citeF-1'></span><span id='citeF-2'></span>[[#cite-1|[1,2]]] and multiphysics simulations <span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-3|[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 <span id='citeF-5'></span><span id='citeF-6'></span>[[#cite-5|[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 <span id='citeF-7'></span><span id='citeF-8'></span>[[#cite-7|[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 <span id='citeF-9'></span><span id='citeF-10'></span>[[#cite-9|[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. <span id='citeF-11'></span>[[#cite-11|[11]]] extends the investigation of <span id='citeF-7'></span>[[#cite-7|[7]]] 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. <span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span>[[#cite-12|[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 <span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-17|[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 <span id='citeF-20'></span>[[#cite-20|[20]]]. 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 <span id='citeF-21'></span>[[#cite-21|[21]]], 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. <span id='citeF-22'></span>[[#cite-22|[22]]]. 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) <span id='citeF-23'></span>[[#cite-23|[23]]]. 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 <span id='citeF-24'></span><span id='citeF-25'></span>[[#cite-24|[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.

==2. Mathematical model==

We consider an arbitrary body, <math display="inline">\Omega \in \mathbb{R}^{\hbox{dim}}</math>, with boundary <math display="inline">\partial \Omega </math>. The displacement of a point <math display="inline">\mathbf{x}\in \Omega </math> is denoted by <math display="inline">\boldsymbol{u}(\mathbf{x}) \in \mathbb{R}^{\hbox{dim}}</math>. The subscript brackets <math display="inline">(\cdot )_{[]}</math> indicates the component of the vector or matrix. By assuming small deformations and deformation gradients, the infinitesimal strain tensor, denoted <math display="inline">\boldsymbol{\varepsilon }(\mathbf{x}) \in \mathbb{R}^{\hbox{dim}\times \hbox{dim}}</math>, is given by

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\varepsilon }= \frac{1}{2}\left(\nabla \boldsymbol{u}+ \left[\nabla \boldsymbol{u}\right]^T \right). </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

Moreover, by considering isotropic elastic materials, the stress tensor, <math display="inline">\boldsymbol{\sigma }(\mathbf{x}) \in \mathbb{R}^{\hbox{dim}\times \hbox{dim}}</math>, is defined as

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }= 2\mu \boldsymbol{\varepsilon }+ \lambda tr(\boldsymbol{\varepsilon }) \mathbf{I}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

where <math display="inline">\mathbf{I}</math> is the identity matrix, <math display="inline">tr(\cdot )</math> is the trace, and <math display="inline">\lambda </math> and <math display="inline">\mu </math> are the Lamé parameters characterizing the material. These parameters are related with the Young's modulus, <math display="inline">E</math>, and the Poisson's ratio, <math display="inline">\nu </math>, by the following equivalences

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu = \frac{E}{2(1+\mu )}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\lambda = \frac{\nu E}{(1+\nu )(1- \lambda _\nu )}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where <math display="inline">\lambda _{\nu } = \nu </math> for plane stress analysis, and <math display="inline">\lambda _{\nu } = 2\nu </math> for plane strain and the 3D cases.

The strong form of the static linear elasticity equation is

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla \cdot \boldsymbol{\sigma }= 0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

The domain boundary <math display="inline">\partial \Omega = \Gamma _N \cup \Gamma _D</math> is the union of the boundary with Neumann conditions, denoted <math display="inline">\Gamma _N</math>, and the boundary with Dirichlet conditions, denoted <math display="inline">\Gamma _D</math>. The equation [[#eq-5|(5)]] must be satisfied for the given forces <math display="inline">\boldsymbol{b}_N</math> and displacements <math display="inline">\boldsymbol{u}_D</math> along the boundary. The following equations remark these user defined boundary conditions,

<span id="eq-6.a"></span>
<span id="eq-6.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }\mathbf{n}= \boldsymbol{b}_N(\mathbf{x}),  \mathbf{x}\in \Gamma _N,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.a)
|-
| style="text-align: center;" | <math>  \boldsymbol{u}= \boldsymbol{u}_D(\mathbf{x}),   \mathbf{x}\in \Gamma _D, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (6.b)
|}
|}

The Figure [[#img-1|1]] illustrates the initial body with the boundary conditions, and the distorted body after equilibrium is solved in the elasticity equation.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_4049_image1.jpg|600px|(a) The initial body, Ω, with its boundary conditions   ퟃΩ= Γ<sub>N</sub>∪Γ<sub>D</sub>. (b) The distorted body            resulting from solving equilibrium in the elasticity equation,            with the boundary conditions given by b<sub>N</sub> and u<sub>D</sub>.            The boundaries where there are not conditions indicated explicitly,            correspond to Neumann conditions b<sub>N</sub>= 0.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1'''. (a) The initial body, <math>\Omega </math>, with its boundary conditions   <math>\partial \Omega = \Gamma _N \cup \Gamma _D</math>. (b) The distorted body            resulting from solving equilibrium in the elasticity equation,            with the boundary conditions given by <math>\boldsymbol{b}_N</math> and <math>\boldsymbol{u}_D</math>.  The boundaries where there are not conditions indicated explicitly, correspond to Neumann conditions <math>\boldsymbol{b}_N = 0</math>
|}


==3. Numerical method==

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.

===3.1 Discretization of domain into control volumes===

The domain <math display="inline">\Omega </math>  is discretized into <math display="inline">N</math> control volumes, denoted <math display="inline">V_i</math>, using the Control Volume Function Approximation (CVFA) proposed by Li et al. <span id='citeF-24'></span><span id='citeF-25'></span>[[#cite-24|[24,25]]]. The partition <math display="inline">P_h</math> of <math display="inline">\Omega </math> is defined by

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>P_h = \bigcup ^N_{i=1} V_i, ~~~\hbox{with}~~~  V_i\cap V_j = \emptyset , ~~~ i \neq j, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where the boundary of each control volume, <math display="inline">\partial V_i</math>, is composed by <math display="inline">N_i</math> flat faces, denoted <math display="inline">e_{ij}</math>,

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial V_i= \bigcup ^{N_i}_{j=1} e_{ij}, ~~~\hbox{with}~~~  e_{ij}\cap e_{ik} = \emptyset ,~~~j \neq k. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

The Figure [[#img-2|2]] illustrates the partition <math display="inline">P_h</math> of <math display="inline">\Omega </math> into <math display="inline">N</math> control volumes defined in the equations [[#eq-7|(7)]] and [[#eq-8|(8)]]. <div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_8669_image2.jpg|600px|The partition Pₕ is the discretization of the domain Ω into   N control volumes.            The boundary of the control volumes, ퟃV<sub>i</sub>, is            conformed by N<sub>i</sub> flat faces, denoted e<sub>ij</sub>. The unit vector   n<sub>ij</sub> is normal to the face e<sub>ij</sub>.            The faces of the volumes adjacent to the boundary Γ<sub>N</sub> are            integrated using the condition b<sub>N</sub>.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2'''. The partition <math>P_h</math> is the discretization of the domain <math>\Omega </math> into   <math>N</math> control volumes.            The boundary of the control volumes, <math>\partial V_i</math>, is            conformed by <math>N_{i}</math> flat faces, denoted <math>e_{ij}</math>. The unit vector   <math>\mathbf{n}_{ij}</math> is normal to the face <math>e_{ij}</math>.            The faces of the volumes adjacent to the boundary <math>\Gamma _N</math> are            integrated using the condition <math>\boldsymbol{b}_N</math>
|}


The Figure [[#img-3|3]] shows a three dimensional control volume. 

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_1413_image3.jpg|600px|The boundary ퟃV<sub>i</sub> of the three dimensional control volume V<sub>i</sub> is            subdivided into N<sub>i</sub> flat faces, denoted e<sub>ij</sub>. The unit vector   n<sub>ij</sub> is normal to the face e<sub>ij</sub>.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3'''. The boundary <math>\partial V_i</math> of the three dimensional control volume <math>V_i</math> is            subdivided into <math>N_{i}</math> flat faces, denoted <math>e_{ij}</math>. The unit vector   <math>\mathbf{n}_{ij}</math> is normal to the face <math>e_{ij}</math>
|}

Every control volume <math display="inline">V_i</math> must have a calculation point

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}_i \in V_i\cup \partial V_i, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

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 <math display="inline">\Gamma _D</math>, it is convenient to establish the calculation point over the corresponding boundary face,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}_i \in \partial V_i\cap \Gamma _D, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

in order to set the Dirichlet condition directly on the point.

===3.2 Control volumes integration===

The elasticity equation [[#eq-5|(5)]] is now integrated over the control volume

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{V_i} \nabla \cdot \boldsymbol{\sigma }~dV = 0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

using the divergence theorem we transform the volume integral into a surface  integral over the volume boundary

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{\partial V_i} \boldsymbol{\sigma }\mathbf{n}~dS = 0. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

The evaluation of the surface integrals is based on the approximation of the displacement field inside the neighborhood of the volume, denoted <math display="inline">{B}_i</math>,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{u}_i(\mathbf{x}) = \sum \limits _{q\in {B}_i}\varphi _q \mathbf{x}_q, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

making use of a group of piece-wise polynomial interpolators, denoted <math display="inline">\varphi _q</math>. 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 [[#eq-1|(1)]] and stress [[#eq-2|(2)]] definitions. Taking advantage of the stress tensor symmetry <math display="inline">~\boldsymbol{\sigma }</math>, we rewrite the stress normal to the boundary as

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }\mathbf{n}= \begin{bmatrix}\boldsymbol{\sigma }_{[11]} & \boldsymbol{\sigma }_{[12]} \\ \boldsymbol{\sigma }_{[12]} & \boldsymbol{\sigma }_{[22]} \end{bmatrix}  \begin{bmatrix}\mathbf{n}_{[1]} \\ \mathbf{n}_{[2]} \end{bmatrix} = \begin{bmatrix}\mathbf{n}_{[1]} & & \mathbf{n}_{[2]} \\  & \mathbf{n}_{[2]} & \mathbf{n}_{[1]} \end{bmatrix} \begin{bmatrix}\boldsymbol{\sigma }_{[11]} \\ \boldsymbol{\sigma }_{[22]} \\ \boldsymbol{\sigma }_{[12]} \end{bmatrix} = \mathbf{T}\vec{\boldsymbol{\sigma }}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">\mathbf{T}</math> is the face orientation matrix and <math display="inline">~\vec{\boldsymbol{\sigma }}</math> is the engineering stress vector. Developing the stress definition [[#eq-2|(2)]] component-wise we can decompose it into the constitutive matrix, denoted <math display="inline">\mathbf{D}</math>, and the engineering strain vector, denoted <math display="inline">\vec{\boldsymbol{\varepsilon }}</math>, as follows

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vec{\boldsymbol{\sigma }} =  \begin{bmatrix}\boldsymbol{\sigma }_{[11]} \\ \boldsymbol{\sigma }_{[22]} \\ \boldsymbol{\sigma }_{[12]} \end{bmatrix} = \begin{bmatrix}2\mu ~ \boldsymbol{\varepsilon }_{[11]} + \lambda \left(\boldsymbol{\varepsilon }_{[11]} + \boldsymbol{\varepsilon }_{[22]}\right)\\ 2\mu ~ \boldsymbol{\varepsilon }_{[22]} + \lambda \left(\boldsymbol{\varepsilon }_{[11]} + \boldsymbol{\varepsilon }_{[22]}\right)\\ 2\mu ~ \boldsymbol{\varepsilon }_{[12]} \end{bmatrix} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|-
| style="text-align: center;" | <math>  =\begin{bmatrix}\left(2\mu + \lambda \right)& \lambda &  \\ \lambda & \left(2\mu + \lambda \right)& \\ & & \mu  \end{bmatrix} \begin{bmatrix}\boldsymbol{\varepsilon }_{[11]} \\ \boldsymbol{\varepsilon }_{[22]} \\ 2\boldsymbol{\varepsilon }_{[12]} \end{bmatrix} = \mathbf{D}\vec{\boldsymbol{\varepsilon }}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}
|}

then the components of the strain vector are retrieved from the equation [[#eq-1|(1)]], and it is decomposed into the matrix differential operator <math display="inline">\mathbf{S}</math> and the displacement function <math display="inline">\boldsymbol{u}</math>.

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\vec{\boldsymbol{\varepsilon }} = \begin{bmatrix}\boldsymbol{\varepsilon }_{[11]} \\ \boldsymbol{\varepsilon }_{[22]} \\ 2\boldsymbol{\varepsilon }_{[12]} \end{bmatrix} = \begin{bmatrix}\displaystyle \frac{\partial \boldsymbol{u}_{[1]}}{\partial \mathbf{x}_{[1]}}\\ \displaystyle \frac{\partial \boldsymbol{u}_{[2]}}{\partial \mathbf{x}_{[2]}}\\ \displaystyle \frac{\partial \boldsymbol{u}_{[1]}}{\partial \mathbf{x}_{[2]}} +  \displaystyle \frac{\partial \boldsymbol{u}_{[2]}}{\partial \mathbf{x}_{[1]}} \end{bmatrix} = \begin{bmatrix}\displaystyle  \frac{\partial }{\partial \mathbf{x}_{[1]}} &\\ & \displaystyle \frac{\partial }{\partial \mathbf{x}_{[2]}}\\ \displaystyle  \frac{\partial }{\partial \mathbf{x}_{[2]}} & \displaystyle  \frac{\partial }{\partial \mathbf{x}_{[1]}} \end{bmatrix} \begin{bmatrix}\boldsymbol{u}_{[1]} \\ \boldsymbol{u}_{[2]} \end{bmatrix} = \mathbf{S}\boldsymbol{u}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

Summarizing the equations [[#eq-14|(14)]], [[#eq-16|(16)]] and [[#eq-17|(17)]] we have

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }\mathbf{n}= \mathbf{T}\vec{\boldsymbol{\sigma }} = \mathbf{T}\mathbf{D}\vec{\boldsymbol{\varepsilon }} = \mathbf{T}\mathbf{D}\mathbf{S}\boldsymbol{u}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

where <math display="inline">\mathbf{T}\mathbf{D}\mathbf{S}</math> is the stiffness of the volume boundary.

Once the displacement field is decoupled, we rewrite the equation [[#eq-12|(12)]] as

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{\partial V_i} \mathbf{T}\mathbf{D}\mathbf{S}\boldsymbol{u}~dS = 0. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

Using the fact that the control volume boundary is divided into flat faces, as in equation [[#eq-8|(8)]], we split the integral [[#eq-19|(19)]] into the sum of the flat faces integrals

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _{j=1}^{N_i} \int _{e_{ij}} \mathbf{T}\mathbf{D}\mathbf{S}\boldsymbol{u}~dS = 0. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

Notice that the face orientation <math display="inline">\mathbf{T}</math> along the flat face, denoted <math display="inline">\mathbf{T}_{ij}</math>, 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 <math display="inline">\mathbf{D}</math> along the flat face, denoted <math display="inline">\mathbf{D}_{ij}</math>, is also considered constant. The matrix <math display="inline">\mathbf{D}_{ij}</math> is estimated from the harmonic average of the Lamé parameters assigned to the adjacent volumes, where <math display="inline">\lambda _i</math> and <math display="inline">\mu _i</math> are the properties of the volume <math display="inline">V_i</math>,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mu _{ij} = \frac{2\mu _i\mu _j}{\mu _i + \mu _j} ~~~~\hbox{and}~~~~ \lambda _{ij} = \frac{2\lambda _i\lambda _j}{\lambda _i + \lambda _j}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

with <math display="inline">\mathbf{T}_{ij}</math> and <math display="inline">\mathbf{D}_{ij}</math> we simplify the equation [[#eq-20|(20)]] as

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _{j=1}^{N_i} \mathbf{T}_{ij}\mathbf{D}_{ij}\mathbf{H}_{ij}= 0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

where <math display="inline">\mathbf{H}_{ij}</math> is the strain integral along the flat face <math display="inline">e_{ij}</math>,

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{H}_{ij}= \int _{e_{ij}} \mathbf{S}\boldsymbol{u}~dS. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

The accuracy of the method depends on the correct evaluation of this integral.

===3.3 Calculating face integrals===

The surface integrals <math display="inline">\mathbf{H}_{ij}</math> along the flat faces <math display="inline">e_{ij}</math> 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 <math display="inline">\mathbf{x}_i</math> from the neighborhood of <math display="inline">V_i</math>. 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 <math display="inline">{B}_i</math> of volume <math display="inline">V_i</math> as the minimum set of calculation points <math display="inline">\mathbf{x}_j</math> such that the simplices intersecting <math display="inline">V_i</math> does not change if we add another calculation point to the set. Observe that the neighborhood <math display="inline">{B}_i</math> does not always coincide with the set of calculation points in volumes adjacent to <math display="inline">V_i</math>, as in most of the FV formulations. Once the neighborhood <math display="inline">{B}_i</math> is triangulated, we ignore the simplices with angles smaller than <math display="inline">10</math> degrees, and the simplices formed outside the domain, which commonly appear in concavities of the boundary <math display="inline">\partial \Omega </math>. The local set of simplices resulting from the neighborhood of <math display="inline">V_i</math> is denoted <math display="inline">P_\alpha </math>. The Figure [[#img-4|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 <math display="inline">{B}_i</math>.

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_7534_image4.jpg|600px|(a) The dotted line illustrates the triangulation of the calculation            points of adjacent volumes to V<sub>i</sub>, used by most of the FV methods.            (b) The dotted line shows the simplices forming the piece-wise            approximation used to solve the integrals H<sub>ij</sub> of the            control volume V<sub>i</sub>.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4'''. (a) The dotted line illustrates the triangulation of the calculation            points of adjacent volumes to <math>V_i</math>, used by most of the FV methods. (b) The dotted line shows the simplices forming the piece-wise approximation used to solve the integrals <math>\mathbf{H}_{ij}</math> of the control volume <math>V_i</math>
|}


The face <math display="inline">e_{ij}</math> is subdivided into <math display="inline">N_{ij}</math> subfaces, denoted <math display="inline">e_{ijk}</math>,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e_{ij}= \bigcup \limits _{k=1}^{N_{ij}} e_{ijk},   ~~~\hbox{with}~~~   e_{ijk}\cap e_{ijl} = \emptyset ,~~~   l \neq k, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

these subfaces result from the intersection between <math display="inline">P_\alpha </math> and the control volume <math display="inline">V_i</math>. The Figure [[#img-4|4]].b illustrates six key points of this approach, 1) the simplices are used to create a polynomial interpolation of <math display="inline">\boldsymbol{u}(\mathbf{x})</math> 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 <math display="inline">V_i</math> and <math display="inline">V_k</math>, 4) there are volumes that require information of non-adjacent volumes to calculate its face integrals, such as <math display="inline">V_i</math> requires <math display="inline">V_k</math>, 5) the dependency between volumes is not always symmetric, which means that if <math display="inline">V_i</math> requires <math display="inline">V_k</math> does not implies that <math display="inline">V_k</math> requires <math display="inline">V_i</math>, and 6) non conforming meshes are supported, as shown in the faces formed by <math display="inline">V_a</math>, <math display="inline">V_b</math>, <math display="inline">V_c</math>, <math display="inline">V_d</math> and <math display="inline">V_j</math>.

The integral [[#eq-23|(23)]] is now rewritten in terms of the subfaces

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{H}_{ij}= \sum \limits _{k=1}^{N_{ij}}~\int _{e_{ijk}} \mathbf{S}\boldsymbol{u}~dS, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

Each subface <math display="inline">e_{ijk}</math> is bounded by a simplex, where the displacement <math display="inline">\boldsymbol{u}_{ijk}</math>, and it derivatives, <math display="inline">~(\mathbf{S}\boldsymbol{u})_{ijk}</math>, are a polynomial interpolation. Hence the integrals in equation [[#eq-25|(25)]] are solved exactly by using the Gauss-Legendre quadrature with the required number of integration points, denoted <math display="inline">N_g</math>, depending on the polynomial degree,

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{e_{ijk}} \mathbf{S}\boldsymbol{u}~dS = \sum \limits _{g = 1}^{N_{g}} ~~   w_g (\mathbf{S}\boldsymbol{u})_{ijk}|_{\mathbf{x}_g}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

where <math display="inline">w_g</math> is the corresponding quadrature weight and <math display="inline">(\mathbf{S}\boldsymbol{u})_{ijk}|_{\mathbf{x}_g}</math> is the strain evaluation of the Gauss point with the proper change of interval, denoted <math display="inline">\mathbf{x}_{g}</math>. The Figure [[#img-5|5]] shows the change of interval required for a 2D face. A 3D face (a polygon) must be subdivided to be integrated with a triangular quadrature. <div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_3659_image5.jpg|600px|(a) The shaded volume V<sub>i</sub> is being integrated.            The integral over the subface e<sub>ijk</sub> is calculated using the            polynomial approximation of the shaded simplex.            The integration point must be mapped to            (b) the normalized space [-1,1] in order to use the Gauss-            Legendre quadrature.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5'''. (a) The shaded volume <math>V_i</math> is being integrated.            The integral over the subface <math>e_{ijk}</math> is calculated using the            polynomial approximation of the shaded simplex.            The integration point must be mapped to            (b) the normalized space <math>[-1,1]</math> in order to use the Gauss-            Legendre quadrature
|}


Most of the cases, the displacement <math display="inline">~\boldsymbol{u}_{ijk}</math> is interpolated inside the simplices, but in some geometrical locations these can not be created, in consequence, the displacement <math display="inline">~\boldsymbol{u}_{ijk}</math> is interpolated pair-wise using the volumes adjacent to the subface <math display="inline">~e_{ijk}</math>. We discuss both strategies in the following subsections.

===3.4 Simplex-wise polynomial approximation===

In the general case, the simplices are formed by <math display="inline">~(\hbox{dim}+1)~</math> points. The points forming the simplex that is bounding the subface <math display="inline">e_{ijk}</math> are denoted <math display="inline">\mathbf{x}_q</math>, and its displacements <math display="inline">\boldsymbol{u}_q</math>.

The shape functions used for the polynomial interpolation are defined into the normalized space. A point in such space is denoted <math display="inline">\boldsymbol{\xi }</math>, its <math display="inline">d^{th}</math> component is denoted <math display="inline">\boldsymbol{\xi }_{[d]}</math>, and the <math display="inline">q^{th}</math> point forming the simplex is expressed <math display="inline">\boldsymbol{\xi }_q</math>. The nodes of the normalized simplex are given by the origin, <math display="inline">\boldsymbol{0}</math>, and the standard basis vectors,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\xi }_q = \left\{ \begin{array}{rl}\mathbf{e}_q, & \hbox{ for} ~~ q \in [1,\hbox{ dim}],\\ \boldsymbol{0}, & \hbox{ if} ~~ q=\hbox{ dim}+1 \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline">\mathbf{e}_q</math> is the <math display="inline">q^{th}</math> standard basis vector. The Figures [[#img-6|6]] and [[#img-7|7]] illustrates the original and the normalized simplices with the corresponding node numeration for 2D and 3D respectively. <div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_1743_image6.jpg|600px|(a) The simplex formed by the points x₁, x₂ and   x₃ in the original space contains an interior point   x<sub>g</sub> that is mapped to            (b) ξ<sub>g</sub> into the normalized 2D-simplex formed by the points   ξ₁, ξ₂ and ξ₃.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6'''. (a) The simplex formed by the points <math>\mathbf{x}_1</math>, <math>\mathbf{x}_2</math> and   <math>\mathbf{x}_3</math> in the original space contains an interior point   <math>\mathbf{x}_g</math> that is mapped to            (b) <math>\boldsymbol{\xi }_g</math> into the normalized 2D-simplex formed by the points   <math>\boldsymbol{\xi }_1</math>, <math>\boldsymbol{\xi }_2</math> and <math>\boldsymbol{\xi }_3</math>
|}


<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_8671_image7.jpg|600px|(a) The 3D-simplex formed by the points x₁, x₂,   x₃ and x₄  in the original space contains an interior            point x<sub>g</sub> that is mapped to            (b) ξ<sub>g</sub> into the normalized 3D-simplex formed by the points   ξ₁, ξ₂, ξ₃ and ξ₄.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7'''. (a) The 3D-simplex formed by the points <math>\mathbf{x}_1</math>, <math>\mathbf{x}_2</math>,   <math>\mathbf{x}_3</math> and <math>\mathbf{x}_4</math>  in the original space contains an interior            point <math>\mathbf{x}_g</math> that is mapped to            (b) <math>\boldsymbol{\xi }_g</math> into the normalized 3D-simplex formed by the points   <math>\boldsymbol{\xi }_1</math>, <math>\boldsymbol{\xi }_2</math>, <math>\boldsymbol{\xi }_3</math> and <math>\boldsymbol{\xi }_4</math>
|}


The shape functions, denoted <math display="inline">\varphi _q</math>, are used to interpolate the displacement field inside the normalized simplex. Such functions are non-negative and are given by

<span id="eq-28.a"></span>
<span id="eq-28.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{P}_c\left(\boldsymbol{\xi }_{[q]}\right),  ~~~\hbox{ if} ~~ q \in [1,\hbox{ dim}],       </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (28.a)
|-
| style="text-align: center;" | <math>              1 - \sum \limits _{d=1}^{\hbox{dim}}       {P}_c\left(\boldsymbol{\xi }_{[d]}\right),  ~~~\hbox{ for} ~~       q=\hbox{ dim}+1,   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (28.b)
|}
|}

where <math display="inline">{P}_c(\cdot )</math> is the ''homeostatic spline'', which is the simplest polynomial defined in the interval <math display="inline">[0,1]</math> that have <math display="inline">~c~</math> 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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _{q=1}^{\hbox{dim} + 1}   \varphi _q(\boldsymbol{\xi }) = 1 ~~~~   \hbox{for any }\boldsymbol{\xi }\hbox{ inside the simplex,} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

furthermore, these functions are equal to one in its corresponding node, which implies that

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varphi _q(\boldsymbol{\xi }_q)  = 1  ~~~~   \hbox{ for any }\boldsymbol{\xi }_q\hbox{ forming the simplex,}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|-
| style="text-align: center;" | <math>   \varphi _q(\boldsymbol{\xi }_p)  = 0  ~~~~   \hbox{ for any }\boldsymbol{\xi }_p \neq \boldsymbol{\xi }_q\hbox{ forming the simplex,} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}
|}

The gradients of the shape functions with respect to the normalized space are denoted <math display="inline">\nabla _{\boldsymbol{\xi }}\varphi _q</math>. The norm of the sum of such gradients is zero

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left|\left|\sum \limits _{q=1}^{\hbox{dim} + 1}   \nabla _{\boldsymbol{\xi }}\varphi _q\left(\boldsymbol{\xi }\right)   \right|\right|= 0 ~~~~   \hbox{for any }\boldsymbol{\xi }\hbox{ inside the simplex,} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

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, <math display="inline">\boldsymbol{p}\left(\boldsymbol{\xi }\right)</math>, by using the shape functions and the points forming the simplex

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{p}\left(\boldsymbol{\xi }\right)=   \sum \limits _{q=1}^{\hbox{dim}+1} \varphi _q\left(\boldsymbol{\xi }\right)   \mathbf{x}_q, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

In order to calculate the normalized point, denoted <math display="inline">\boldsymbol{\xi }_g</math>, associated to the integration point <math display="inline">\mathbf{x}_g = \boldsymbol{p}\left(\boldsymbol{\xi }_g\right)</math>, we use the shape functions definitions to rewrite the equation [[#eq-33|(33)]] in matrix form

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{p}\left(\boldsymbol{\xi }\right)=   \underbrace{   \begin{bmatrix}     \mathbf{x}_{3[1]} \\ \mathbf{x}_{3[2]}   \end{bmatrix} +   \left[     \begin{matrix}       \left(\mathbf{x}_{1[1]} - \mathbf{x}_{3[1]}\right)&       \left(\mathbf{x}_{2[1]} - \mathbf{x}_{3[1]}\right)\\       \left(\mathbf{x}_{1[2]} - \mathbf{x}_{3[2]}\right)&       \left(\mathbf{x}_{2[2]} - \mathbf{x}_{3[2]}\right)     \end{matrix}   \right]   \begin{bmatrix}     {P}_c\left(\boldsymbol{\xi }_{[1]}\right)\\ {P}_c\left(\boldsymbol{\xi }_{[2]}\right)   \end{bmatrix}   }_{\hbox{2D case (triangle)}}</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|-
| style="text-align: center;" | <math>      =\mathbf{x}_{(\hbox{dim}+1)} + \mathbf{J}_{\Delta }~{P}_c\left(\boldsymbol{\xi }\right), </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}
|}

where <math display="inline">{P}_c\left(\boldsymbol{\xi }\right)</math> is the vector resulting from evaluating the spline for <math display="inline">\boldsymbol{\xi }</math> component-wise, and <math display="inline">\mathbf{J}_{\Delta }~</math> is the distortion matrix given by the concatenation of the following column vector differences

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{J}_{\Delta }~= \left[     (\mathbf{x}_1 - \mathbf{x}_{(\hbox{dim}+1)}), ... ,      (\mathbf{x}_{\hbox{dim}} - \mathbf{x}_{(\hbox{dim}+1)})     \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Now, from equation [[#eq-35|(35)]] we retrieve the point <math display="inline">\mathbf{x}_g</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}_g = \boldsymbol{p}\left(\boldsymbol{\xi }_g\right)=   \mathbf{x}_{(\hbox{dim}+1)} + \mathbf{J}_{\Delta }~{P}_c\left(\boldsymbol{\xi }_g\right), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

and solving for <math display="inline">\boldsymbol{\xi }_g</math> we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\xi }_g = {Q}_c   \left(   \left(\mathbf{J}_{\Delta }~\right)^{-1}   \left(\mathbf{x}_g - \mathbf{x}_{(\hbox{dim}+1)} \right)   \right), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

where <math display="inline">{Q}_c</math> is the inverse function of <math display="inline">{P}_c</math> applied component-wise to the product of the matrix-vector operation.

Similar to the approximation in equation [[#eq-33|(33)]], within the simplex enclosing the subface <math display="inline">e_{ijk}</math>, the displacement field evaluated at <math display="inline">\mathbf{x}_g</math> is defined as,

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{u}_{ijk}|_{\mathbf{x}_g} =   \sum \limits _{q=1}^{\hbox{dim} + 1}   \varphi _{q}\left(\boldsymbol{\xi }_g\right)~\boldsymbol{u}_q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

Hence, when calculating the quadrature of equation [[#eq-26|(26)]], the strain evaluated at the integration point is given by

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\mathbf{S}\boldsymbol{u}\right)_{ijk}|_{\mathbf{x}_g}~~~ =   \sum \limits _{q=1}^{\hbox{dim} + 1}   \mathbf{S}\varphi _q(\boldsymbol{\xi }_g)   ~~\boldsymbol{u}_q,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|-
| style="text-align: center;" | <math>   =\left[   \begin{matrix}\displaystyle \frac{\partial \varphi _1}{\partial \mathbf{x}_{[1]}} & \\     & \displaystyle \frac{\partial \varphi _1}{\partial \mathbf{x}_{[2]}} \\     \displaystyle \frac{\partial \varphi _1}{\partial \mathbf{x}_{[2]}} &     \displaystyle \frac{\partial \varphi _1}{\partial \mathbf{x}_{[1]}}   \end{matrix}~~~~~   \begin{matrix}\displaystyle \frac{\partial \varphi _2}{\partial \mathbf{x}_{[1]}} & \\     & \displaystyle \frac{\partial \varphi _2}{\partial \mathbf{x}_{[2]}} \\     \displaystyle \frac{\partial \varphi _2}{\partial \mathbf{x}_{[2]}} &     \displaystyle \frac{\partial \varphi _2}{\partial \mathbf{x}_{[1]}}   \end{matrix}~~~~~   \begin{matrix}\displaystyle \frac{\partial \varphi _3}{\partial \mathbf{x}_{[1]}} & \\     & \displaystyle \frac{\partial \varphi _3}{\partial \mathbf{x}_{[2]}} \\     \displaystyle \frac{\partial \varphi _3}{\partial \mathbf{x}_{[2]}} &     \displaystyle \frac{\partial \varphi _3}{\partial \mathbf{x}_{[1]}}   \end{matrix}   \right]_{|\mathbf{x}_g}   \begin{bmatrix}\boldsymbol{u}_{1[1]} \\ \boldsymbol{u}_{1[2]} \\     \boldsymbol{u}_{2[1]} \\ \boldsymbol{u}_{2[2]} \\     \boldsymbol{u}_{3[1]} \\ \boldsymbol{u}_{3[2]}   \end{bmatrix} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|-
| style="text-align: center;" | 
|-
| style="text-align: center;" | <math>   =\mathbf{B}_{ijk}|_{\mathbf{x}_g} ~~~ \vec{\boldsymbol{u}}_{ijk}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}
|}

where <math display="inline">\mathbf{B}_{ijk}|_{\mathbf{x}_g}</math> captures the deformation at <math display="inline">\mathbf{x}_g</math>, and <math display="inline">\vec{\boldsymbol{u}}_{ijk}</math> is the vector with the concatenated displacement components of the points forming the simplex.

In order to calculate the deformation matrix <math display="inline">\mathbf{B}_{ijk}</math>, we require the derivatives of the shape functions with respect to <math display="inline">\mathbf{x}</math>, denoted <math display="inline">\nabla \varphi _q</math>. These derivatives are calculated by solving the linear systems resulting from the chain rule

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla _{\boldsymbol{\xi }}\varphi _{q} =  \begin{bmatrix}\displaystyle \frac{\partial \varphi _q}{\partial \boldsymbol{\xi }_{[1]}} \\ \displaystyle \frac{\partial \varphi _q}{\partial \boldsymbol{\xi }_{[2]}} \end{bmatrix} = \begin{bmatrix}\displaystyle \frac{\partial \mathbf{x}_{[1]}}{\partial \boldsymbol{\xi }_{[1]}} & \displaystyle \frac{\partial \mathbf{x}_{[2]}}{\partial \boldsymbol{\xi }_{[1]}} \\ \displaystyle \frac{\partial \mathbf{x}_{[1]}}{\partial \boldsymbol{\xi }_{[2]}} & \displaystyle \frac{\partial \mathbf{x}_{[2]}}{\partial \boldsymbol{\xi }_{[2]}} \end{bmatrix} \begin{bmatrix}\displaystyle \frac{\partial \varphi _q}{\partial \mathbf{x}_{[1]}} \\ \displaystyle \frac{\partial \varphi _q}{\partial \mathbf{x}_{[2]}} \end{bmatrix} = \left(\nabla _{\boldsymbol{\xi }}\boldsymbol{p}\right)^T ~\nabla \varphi _q, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

where <math display="inline">\nabla _{\boldsymbol{\xi }}\boldsymbol{p}</math> is the geometric jacobian evaluated at <math display="inline">\boldsymbol{\xi }</math>.  This jacobian relates both spaces, captures the distortion of the simplex, and is derivated from equation [[#eq-33|(33)]],

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\nabla _{\boldsymbol{\xi }}\boldsymbol{p}=   \sum \limits _{q=1}^{\hbox{dim}+1} \mathbf{x}_q \left(\nabla _{\boldsymbol{\xi }}\varphi _q\right)^T, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

The gradients of the shape functions with respect to <math display="inline">\boldsymbol{\xi }</math> inside the sum are obtained straightforward once we have the spline first derivative <math display="inline">{P}_c'</math>. Notice that the geometric jacobian is equivalent to the distortion matrix <math display="inline">\mathbf{J}_{\Delta }~</math> if and only if the homeostatic spline is <math display="inline">{P}_c(z) = z</math>.

===3.5 Pair-wise polynomial approximation===

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. Figure [[#img-8|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.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;"
|-
|[[Image:Review_101239056093_9792_image8.jpg|600px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''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 <math display="inline">e_{ijk}</math>  is a pair-wise polynomial approximation between the adjacent volumes, <math display="inline">\mathbf{x}_i</math> and <math display="inline">\mathbf{x}_j</math>, regardless the dimension

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{u}_{ijk}(\mathbf{x}_g) =   \underbrace{     \left(1 - {P}_c\left(z_{g}\right)\right)     }_{\varphi _i}  ~\boldsymbol{u}_i +   \underbrace{{P}_c\left(z_{g}\right)}_{\varphi _j} ~\boldsymbol{u}_j, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

where <math display="inline">\varphi _i</math> and <math display="inline">\varphi _j</math> are the shape functions, and <math display="inline">z_{g}</math> is the normalized projection of the integration point <math display="inline">\mathbf{x}_g</math> into the vector which goes from <math display="inline">\mathbf{x}_i</math> to <math display="inline">\mathbf{x}_j</math>, denoted <math display="inline">\mathbf{x}_{\vec{ij}} = \left(\mathbf{x}_j - \mathbf{x}_i\right),</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>z_{g} = \frac{\left(\mathbf{x}_g - \mathbf{x}_i\right)^T\mathbf{x}_{\vec{ij}}}{||\mathbf{x}_{\vec{ij}}||^2}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

When calculating the quadrature of equation [[#eq-26|(26)]], the pairwise  strain is given by

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\mathbf{S}\boldsymbol{u}\right)_{ijk}|_{\mathbf{x}_g}~~~ =   \mathbf{S}\varphi _i(\boldsymbol{\xi }_g)~\boldsymbol{u}_i ~~+~~   \mathbf{S}\varphi _j(\boldsymbol{\xi }_g)~\boldsymbol{u}_j,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|-
| style="text-align: center;" | <math>   =\left[   \begin{matrix}\displaystyle \frac{\partial \varphi _i}{\partial \mathbf{x}_{[1]}} & \\     & \displaystyle \frac{\partial \varphi _i}{\partial \mathbf{x}_{[2]}} \\     \displaystyle \frac{\partial \varphi _i}{\partial \mathbf{x}_{[2]}} &     \displaystyle \frac{\partial \varphi _i}{\partial \mathbf{x}_{[1]}}   \end{matrix}~~~~~   \begin{matrix}\displaystyle \frac{\partial \varphi _j}{\partial \mathbf{x}_{[1]}} & \\     & \displaystyle \frac{\partial \varphi _j}{\partial \mathbf{x}_{[2]}} \\     \displaystyle \frac{\partial \varphi _j}{\partial \mathbf{x}_{[2]}} &     \displaystyle \frac{\partial \varphi _j}{\partial \mathbf{x}_{[1]}}   \end{matrix}   \right]_{|\mathbf{x}_g}   \begin{bmatrix}\boldsymbol{u}_{i[1]} \\ \boldsymbol{u}_{i[2]} \\     \boldsymbol{u}_{j[1]} \\ \boldsymbol{u}_{j[2]}   \end{bmatrix} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|-
| style="text-align: center;" | 
|-
| style="text-align: center;" | <math>   =\mathbf{B}_{ijk}|_{\mathbf{x}_g} ~~~ \vec{\boldsymbol{u}}_{ijk}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}
|}

In the general case, the gradient is not constant along the face <math display="inline">e_{ij}</math>, since its normal is not necessary aligned with <math display="inline">\mathbf{x}_{\vec{ij}}</math>, as illustrated in Figure [[#img-9|9]]. <div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 70%;"
|-
|[[Image:Review_101239056093_4058_image9.jpg|600px|The gradient of the pairwise approximation is not constant along the            face e<sub>ij</sub>, since its normal is not necessary aligned with   x<sub>\vecij</sub>.            The integration point is projected into x<sub>\vecij</sub> to evaluate            the deformation matrix.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figure 9'''. The gradient of the pairwise approximation is not constant along the            face <math>e_{ij}</math>, since its normal is not necessary aligned with   <math>\mathbf{x}_{\vec{ij}}</math>.            The integration point is projected into <math>\mathbf{x}_{\vec{ij}}</math> to evaluate            the deformation matrix
|}


This pairwise approximation must be used only when necessary because it can not capture the deformation orthogonal to <math display="inline">\mathbf{x}_{\vec{ij}}</math>.

===3.6 Homeostatic spline===

The homeostatic spline is a function of a single variable defined from <math display="inline">z=0</math> to <math display="inline">z=1</math>, denoted <math display="inline">{P}_c(z)</math>, and curved by the parameter <math display="inline">~c</math>, which indicates the level of smoothness. This spline is the simplest polynomial with <math display="inline">~c~</math> derivatives equal to zero at the endpoints <math display="inline">z=0</math> and <math display="inline">z=1</math>. The polynomial degree is given by <math display="inline">2c+1</math>, and such a polynomial requires <math display="inline">N_g = c+1</math> integration points to calculate the exact integral in equation [[#eq-26|(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 <math display="inline">2c+2</math> coefficients of the polynomial. The equations of this system were obtained by forcing the <math display="inline">~c</math> derivatives to be zero at the endpoints. Once we solved the coefficients for the first twenty polynomials, from <math display="inline">c=0</math> to <math display="inline">c=19</math>, we found out that the first half of such coefficients are null, and the entire polynomial can be calculated directly as

<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{P}_c(z) = \frac{1}{b_c}\sum \limits _{k=1}^{1+c} (-1)^{k}~ b_k~~ z^{(k+c)}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

where <math display="inline">b_k</math> is the <math display="inline">k^{th}</math> not null coefficient

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_k = \frac{1}{k+c}   \left(   \prod \limits _{l=1}^{C_k}\frac{\left(1 + c\right)}{l} - 1   \right), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

<math display="inline">C_k</math> is the number of factors needed to calculate <math display="inline">b_k</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>C_k = (c/2) - \big|1 + (c/2) - k\big|, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

and <math display="inline">b_c</math> is accumulation of the coefficients for normalizing the spline

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>b_c = \sum \limits _{k=1}^{1+c}(-1)^{k} b_k, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

The first derivative is simply calculated as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{P}_c'(z) = \frac{1}{b_c}\sum \limits _{k=1}^{1+c} (-1)^{k}~ b_k~~   (k+c) ~~z^{(k+c-1)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

Figure [[#img-10|10]] (a) shows  the evolution of the spline as we increase the smoothness parameter from <math display="inline">c=0</math> to <math display="inline">c=6</math>, and Figure [[#img-10|10]] (b) the evolution of it first derivative. Smoother splines produces higher order polynomials which increases the accuracy of the stress field approximation. This feature is specially important when solving non-linear problems sensibles to the stress field. 

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 65%;"
|-
|[[Image:Review_101239056093_6248_image10.jpg|600px|(a) The evolution of the homeostatic spline from c=0 to c=6            illustrates the smoothness requirements at the endpoints of each            spline and its            (b) first derivatives.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;| '''Figure 10'''. (a) The evolution of the homeostatic spline from <math>c=0</math> to <math>c=6</math> illustrates the smoothness requirements at the endpoints of each spline and its (b) first derivatives
|}


Since the derivatives of the homeostatic spline [[#eq-50|50]] are zero at the endpoints of the interval <math display="inline">[0,1]</math>, the inverse function is not defined in that points. However, we estimate a pseudo-inverse within this interval, <math display="inline">{Q}_c\approx {P}_c^{-1}</math>, by finding the coefficients of a polynomial of the same degree, <math display="inline">2c+1</math>, 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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{Q}_c(0) = {P}_c(0) = 0,~~~ {Q}_c(1) = {P}_c(1) = 1, ~~~\hbox{and}~~~ {Q}_c'(0.5) = \frac{1}{{P}_c'(0.5)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

The higher derivatives in the midpoint are forced to be zero. Once we calculated the coefficients for the first twenty polynomials, from <math display="inline">c=0</math> to <math display="inline">c=19</math>, we found out that the pseudo-inverse can be approximated directly from the following formulae

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{Q}(z) = ~a_1 ~z ~~+~~ (a_1-1)(2c+1)~\sum \limits _{k=1}^{2c}   (-1)^k~ a_k~z^{(k+1)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

where <math display="inline">a_1</math> is the coefficient for <math display="inline">z</math>, and <math display="inline">a_k</math> is the factor that distinguish higher order coefficients. Such terms are calculated as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a_1 = \left(\frac{c}{2\sqrt{2}}+1\right)^2,~~~\hbox{and}~~~ a_k = 2^{(k-1)}~ \prod _{l=1}^{k-1} \left(\frac{2c - l}{2 + l} \right), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

respectively. Figure [[#img-11|11]] exhibits the curves for the first seven levels of smoothness. The null higher derivatives requirement is noticeable at the midpoint.

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 60%;"
|-
|[[Image:Review_101239056093_8102_image11.jpg|700px|Curves of the pseudo-inverse Q<sub>c</sub> for the first seven levels            of smoothness. The slope at the midpoint exposes the null higher            derivatives requirement when increasing the polynomial order.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;| '''Figure 11'''. Curves of the pseudo-inverse <math>{Q}_c</math> for the first seven levels            of smoothness. The slope at the midpoint exposes the null higher            derivatives requirement when increasing the polynomial order
|}


Figure [[#img-12|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 <math display="inline">c=0</math> are the only case where all the shape functions are indistinguishable, these are planes. 

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;"
|-
|[[Image:Review_101239056093_5599_image12.jpg|600px|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.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;| '''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
|}


Figure [[#img-13|13]] shows the magnitude of the gradient with respect to the normalized space. With the same tabular configuration of Figure [[#img-12|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 <math display="inline">~c > 0~</math> the continuity on the stress field is only guaranteed at the calculation points, but not in the simplices edges. 

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;"
|-
|[[Image:Review_101239056093_6320_image13.jpg|600px|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.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;| '''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
|}

===3.7 Assembling volume's equation===

By using the simplex-wise [[#eq-42|(42)]] or the pair-wise [[#eq-49|(49)]] approximation, the strain face integral [[#eq-25|(25)]] is reformulated as

<span id="eq-58"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{H}_{ij}= \sum \limits _{k=1}^{N_{ij}}   \sum \limits _{g=1}^{N_{g}}  ~~w_g ~\mathbf{B}_{ijk}|_{\mathbf{x}_g} ~\vec{\boldsymbol{u}}_{ijk}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

then, the volume equilibrium equation [[#eq-22|(22)]] is

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _{j=1}^{N_i} \mathbf{T}_{ij}\mathbf{D}_{ij}   \sum \limits _{k=1}^{N_{ij}}    \sum \limits _{g=1}^{N_{g}}  ~~w_g ~\mathbf{B}_{ijk}|_{\mathbf{x}_g}   ~\vec{\boldsymbol{u}}_{ijk} = 0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

reordering the terms we obtain

<span id="eq-60"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sum \limits _{j=1}^{N_i}   \sum \limits _{k=1}^{N_{ij}}    \sum \limits _{g=1}^{N_{g}}  ~~w_g ~\mathbf{K}_{ijk}|_{\mathbf{x}_g}   ~\vec{\boldsymbol{u}}_{ijk} = 0, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|}

where the matrix

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{K}_{ijk}|_{\mathbf{x}_g} = \mathbf{T}_{ij}\mathbf{D}_{ij}\mathbf{B}_{ijk}|_{\mathbf{x}_g}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

is the stiffness contribution at the integration point <math display="inline">\mathbf{x}_g</math>, within the subface <math display="inline">e_{ijk}</math>  when integrating the <math display="inline">i^{th}</math> volume. Observe that the stiffness matrix <math display="inline">\mathbf{K}_{ijk}</math> is rectangular and the resulting global stiffness matrix is asymmetric.

===3.8 Boundary conditions===

The Neumann boundary conditions are imposed over the volume faces <math display="inline">e_{ij}</math> intersecting the boundary, by replacing the corresponding term in the sum of equation [[#eq-20|(20)]] with the integral of the function provided in [[#eq-6.a|(6.a)]],

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{e_{ij}} \mathbf{T}\mathbf{D}\mathbf{S}\boldsymbol{u}~dS = \int _{e_{ij}} \boldsymbol{b}_N(\mathbf{x}) ~dS </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

The Dirichlet conditions are imposed over the volumes calculation points by fixing the displacement as it is evaluated in the function given in  [[#eq-6.b|(6.b)]],

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{u}_i = \boldsymbol{u}_D(\mathbf{x}_i), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
|}

Since the Dirichlet conditions are imposed directly on the calculation points, these points must be located along the face <math display="inline">e_{ij}</math> which intersects the  boundary with the condition <math display="inline">\Gamma _D</math>.

===3.9 Special cases===

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 <math display="inline">\mathbf{x}_i</math>. Hence, the subdivision of the neighborhood <math display="inline">{B}_i</math> is already given by the Delaunay triangulation which is dual to the Voronoi mesh, as illustrated in the Figure [[#img-14|14]].a. Moreover, the integrals of subfaces <math display="inline">e_{ijk}</math> using pair-wise approximations can be exactly integrated with the midpoint rule, since the faces are orthogonal to the vector joining the calculation points <math display="inline">\mathbf{x}_{\vec{ij}}</math>, and the derivatives along the subface are constants.

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_3668_image14.jpg|600px|(a) The initial mesh is equivalent to the Voronoi diagram and the            Voronoi centres correspond to the calculation points x<sub>i</sub>.            (b) The initial mesh is generated from a FEM-like triangular mesh.            The calculation points x<sub>i</sub> 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.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14'''. (a) The initial mesh is equivalent to the Voronoi diagram and the            Voronoi centres correspond to the calculation points <math>\mathbf{x}_i</math>.            (b) The initial mesh is generated from a FEM-like triangular mesh.            The calculation points <math>\mathbf{x}_i</math> 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 <math display="inline">\mathbf{x}_i</math> 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 [[#img-14|14]].b. This particular version is equivalent to the cell-centred finite volume scheme introduced by Oñate et al. <span id='citeF-7'></span>[[#cite-7|[7]]], 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.

==4. Results==

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 <span id='citeF-26'></span>[[#cite-26|[26]]]. In such a experiment, the plate is stretched along the horizontal axis with a uniform tension of <math display="inline">\mathbf{f}_{[1]}=10 ~~kPa</math> from each side, as is shown in Figure [[#img-15|15]]. The material is characterized by the Poisson ratio, <math display="inline">\nu=0.3</math>, and Young modulus, <math display="inline">E=10 ~~MPa</math>. Plane stress is assumed with thickness equivalent to the unity. The dimensions of the computational domain are <math display="inline">a=0.5m</math> and <math display="inline">b=2m</math>. 

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_7940_image15.jpg|600px|(a) Infinite plate with a hole being stretched along the horizontal            axis with a force of f<sub>[1]</sub>=10 ~~kPa from each side.            (b) Computational domain, a= 0.5m and b=2m, with axysymmetrical            assumptions used to test the numerical method.            The polar coordinates, r and θ, for calculating the            analytical stress field.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15'''. (a) Infinite plate with a hole being stretched along the horizontal            axis with a force of <math>\mathbf{f}_{[1]}=10 ~~kPa</math> from each side.            (b) Computational domain, <math>a= 0.5m</math> and <math>b=2m</math>, with axysymmetrical            assumptions used to test the numerical method.            The polar coordinates, <math>r</math> and <math>\theta </math>, for calculating the            analytical stress field
|}


The analytical solution is given by the following formulae

<span id="eq-64"></span>
<span id="eq-65"></span>
<span id="eq-66"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }_{[11]} = \mathbf{f}_{[1]}   \left[     1 - \frac{a^2}{r^2}     \left(\frac{3}{2}\cos (2\theta ) + \cos (4\theta )  \right)     + \frac{3a^4}{2 r^4} \cos (4\theta )       \right],</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|-
| style="text-align: center;" | <math>      \boldsymbol{\sigma }_{[22]} = \mathbf{f}_{[1]}   \left[     - \frac{a^2}{r^2}     \left(\frac{1}{2}\cos (2\theta ) - \cos (4\theta )  \right)     - \frac{3a^4}{2 r^4} \cos (4\theta )       \right],</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|-
| style="text-align: center;" | <math>      \boldsymbol{\sigma }_{[12]} = \mathbf{f}_{[1]}   \left[     - \frac{a^2}{r^2}     \left(\frac{1}{2}\sin (2\theta ) + \sin (4\theta )  \right)     + \frac{3a^4}{2 r^4} \sin (4\theta )       \right], </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}
|}

where the polar coordinates,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r = \sqrt{\mathbf{x}_{[1]}^2 + \mathbf{x}_{[2]}^2}, ~~~\hbox{and}~~~ \theta = \tan ^{-1}\left(\frac{\mathbf{x}_{[2]}}{\mathbf{x}_{[1]}}\right), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

are used within the calculus. The Figure [[#img-16|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 <math display="inline">\boldsymbol{\sigma }_{[11]}</math> between <math display="inline">0</math> to <math display="inline">30 ~~kPa</math>, with steps of <math display="inline">1 ~~kPa</math>. (c) Level sets of <math display="inline">\boldsymbol{\sigma }_{[22]}</math> between <math display="inline">-10</math> and <math display="inline">6 ~~kPa</math>, with steps of <math display="inline">0.8 ~~kPa</math>. (d) Level sets of <math display="inline">\boldsymbol{\sigma }_{[12]}</math> between <math display="inline">-10</math> and <math display="inline">2 ~~kPa</math>, with steps of <math display="inline">0.6 ~~kPa</math>. 

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_1864_image16.jpg|600px|(a) Polygonal mesh used for comparison of numerical results.            (b) Level sets of σ<sub>[11]</sub> between 0 to 30 ~~kPa.            (c) Level sets of σ<sub>[22]</sub> between -10 and 6 ~~kPa.            (d) Level sets of σ<sub>[12]</sub> between -10 and 2 ~~kPa.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16'''. (a) Polygonal mesh used for comparison of numerical results.            (b) Level sets of <math>\boldsymbol{\sigma }_{[11]}</math> between <math>0</math> to <math>30 ~~kPa</math>.            (c) Level sets of <math>\boldsymbol{\sigma }_{[22]}</math> between <math>-10</math> and <math>6 ~~kPa</math>.            (d) Level sets of <math>\boldsymbol{\sigma }_{[12]}</math> between <math>-10</math> and <math>2 ~~kPa</math>
|}


The Dirichlet conditions are imposed on the bottom and left side of the computational domain as is shown in the Figure [[#img-16|16]].b. Next in order, the analytic stress of equations [[#eq-64|(64)]], [[#eq-65|(65)]] and [[#eq-66|(66)]] is imposed as Neumann condition over the top and right side of the computational domain.

The Figure [[#img-17|17]].a presents the averaged error as a function of the volumes face length mean, denoted <math display="inline">\Delta x</math>, as we might expect, the error is proportional to the mesh refinement. For a mesh of 478 volumes, the Figure [[#img-17|17]].b shows the percentage of the error with respect to the error of <math display="inline">c=0</math>, for different smoothing levels, <math display="inline">c=0</math> 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. 

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Review_101239056093_1349_image17.jpg|600px|(a) The averaged error decreasing as a function of mesh size, denoted ∆x. (b) Using a mesh of 628 volumes, percentage of error for different smoothing levels with respect to the error of c=0, which is the linear interpolator, error increases after c=2.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17'''. (a) The averaged error decreasing as a function of mesh size, denoted <math>\Delta \mathbf{x}</math>. (b) Using a mesh of 628 volumes, percentage of error for different smoothing levels with respect to the error of <math>c=0</math>, which is the linear interpolator, error increases after <math>c=2</math>
|}

==5. Conclusions==

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.

==Acknowledgement==

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.

==References==

<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

<div id="cite-1"></div>
[[#citeF-1|[1]]] 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.

<div id="cite-2"></div>
[[#citeF-2|[2]]] 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.

<div id="cite-3"></div>
[[#citeF-3|[3]]] 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.

<div id="cite-4"></div>
[[#citeF-4|[4]]] Souli M., Ouahsine A., Lewin L.  ALE formulation for fluid-structure interaction problems. Computer Methods in Applied Mechanics and Engineering, 190:659-675, 2000.

<div id="cite-5"></div>
[[#citeF-5|[5]]] 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.

<div id="cite-6"></div>
[[#citeF-6|[6]]] 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.

<div id="cite-7"></div>
[[#citeF-7|[7]]] 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.

<div id="cite-8"></div>
[[#citeF-8|[8]]] 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.

<div id="cite-9"></div>
[[#citeF-9|[9]]] 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.

<div id="cite-10"></div>
[[#citeF-10|[10]]] 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.

<div id="cite-11"></div>
[[#citeF-11|[11]]] Slone A.K., Bailey C., Cross M. Dynamic solid mechanics using finite volume methods. Applied Mathematical Modeling, 27(2):69-87, 2003.

<div id="cite-12"></div>
[[#citeF-12|[12]]] Demirdzic I., Martinovic D. Finite volume method for thermo-elasto-plastic stress analysis. Computer Methods in Applied Mechanics and Engineering, 109:331-349, 1993.

<div id="cite-13"></div>
[[#citeF-13|[13]]] Demirdzic I., Muzaferija S. Finite volume methods for stress analysis in complex domains. International Journal for Numerical Methods in Engineering, 137:3751-3766, 1994.

<div id="cite-14"></div>
[[#citeF-14|[14]]] 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.

<div id="cite-15"></div>
[[#citeF-15|[15]]] 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.

<div id="cite-16"></div>
[[#citeF-16|[16]]] 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.

<div id="cite-17"></div>
[[#citeF-17|[17]]] 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.

<div id="cite-18"></div>
[[#citeF-18|[18]]] 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.

<div id="cite-19"></div>
[[#citeF-19|[19]]] 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.

<div id="cite-20"></div>
[[#citeF-20|[20]]] 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.

<div id="cite-21"></div>
[[#citeF-21|[21]]] 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.

<div id="cite-22"></div>
[[#citeF-22|[22]]] 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.

<div id="cite-23"></div>
[[#citeF-23|[23]]] Nordbotten J.M. Cell-centered finite volume discretizations for deformable porous media. International Journal for Numerical Methods in Engineering, 100:399-418, 2014.

<div id="cite-24"></div>
[[#citeF-24|[24]]] 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.

<div id="cite-25"></div>
[[#citeF-25|[25]]] 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.

<div id="cite-26"></div>
[[#citeF-26|[26]]] Timoshenko S.P., Goodier J.N. Theory of Elasticity. McGraw-Hill, 3 Edition, 1970.

</div>