Multiscale computational analysis in mechanics using finite calculus. An introduction

<!-- metadata commented in wiki content
==MULTISCALE COMPUTATIONAL ANALYSIS IN MECHANICS USING FINITE CALCULUS. AN INTRODUCTION==

'''Eugenio Oñate

International Center for Numerical Methods in Engineering (CIMNE)

Universidad Politécnica de Cataluña

Edificio C1, Gran Capitán s/n, 08034, Barcelona, Spain

e-mail: onate@cimne.upc.es'''
-->

==Abstract==

The paper introduces a general procedure for computational analysis of a wide class of multiscale problems in mechanics using a finite calculus (FIC) formulation. The FIC approach is based in expressing the governing equations in mechanics accepting that the domain where the standard balance laws are established has a finite size. This introduces naturally additional terms into the classical equations of infinitesimal theory in mechanics which are useful for the numerical solution of problems involving different scales in the physical parameters. The discrete nodal values obtained with the FIC formulation and the finite element method (FEM) can be effectively used as the starting point for obtaining a more refined solution in zones where high gradients of the relevant variables occur using hierarchical or enriched FEM. Typical multiscale problems in mechanics which can be solved with the FIC method include convection-diffusion-reaction problems with high localized gradients, incompressible problems in solid and fluid mechanics, localization problems such as prediction of shear bands in solids and shock waves in compressible fluids, turbulence, etc. The paper presents an introduction of the treatment of multiscale problems using the FIC approach in conjunction with the finite element method (FEM). Examples of application of the FIC/FEM formulation to the solution of simple  multiscale convection-diffusion problems are given.

==1 INTRODUCTION==

Multiscale problems in mechanics are characterized by the existence of numerical values of the geometrical and/or physical parameters differing in several orders of magnitude. Typical examples are convection-diffusion-absorption problems where high values of the convection or the absorption terms lead to sharp boundary and/or internal layers along which the numerical solution can change in many orders of magnitude. A similar problem is found in the analysis of thin discontinuity layers in solids, such as in the case of fracture in concrete, rock or geomaterials, or in shock waves in compressible fluids. A traditional multiscale problem is that of a turbulence fluid, where high gradients of the velocity field occur along random directions of the flow. The transition from a compressible to an incompressible solid or fluid can also be considered as a multiscale problem, in the sense that the propagation speed of sound changes from a finite value (for the compressible case) to an infinite value (in the case of incompressibility). Indeed, problems where many different scales of the material properties co-exist in a solid (such as in composites) are also classical examples of multiscale situations in solid mechanics. For the sake of conciseness, however, the case of composites will not be dealt with in the paper, although it can be also tackled with the general FIC formulation here described.

Recent attempts to develop a unified relativity theory in physics adequate for dealing with the whole spectrum of scales in the universe, from the Plank atomistic scale (<math display="inline">10^{-33}</math> cm) to the cosmological scale (<math display="inline">10^{+28}</math> cm) indicate that a suitable mathematical model should incorporate the effect of the different scales within the governing equations of the model <span id='citeF-1'></span>[[#cite-1|[1]]]. Failing to include these scale terms may lead to unstable behaviour of the equations when solved numerically for very different values of the physical, geometrical and/or time parameters of the  problem.

The Finite Calculus (FIC) method developed by Oñate and co-workers <span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-7'></span>[[#cite-2|[2,3,4,5,6,7]]] is a consistent procedure to re-formulate the governing equation in mechanics introducing new terms involving characteristic space and time dimensions into the equations. The modified equations are derived by invoking the balance laws in mechanics in a space-time domain of finite size. The new terms introduced by the FIC approach are essential to obtain physical (stable) numerical solutions for all ranges of the parameters governing the physical problem.

This paper presents an extension of the FIC method to deal with multiscale problems in mechanics using the finite element method (FEM) <span id='citeF-10'></span>[[#cite-10|[10]]]. The discrete (nodal) values provided by the numerical FEM solution of the FIC equations are taken as the “correct macroscale” solution for a given discretization. A more refined numerical solution in zones where high gradients of the relevant variables occur can be simply obtained by using a hierarchical FE approximations or by introducing enriched interpolation functions accounting for the “correct” shape of the solution sought. The paper describes the basic ideas of the multiscale FIC procedure proposed and two examples of application to simple 1D and 2D convection-diffusion problems showing the applicability of the method. The possibilities of the FIC method for strain localization problems are also briefly highlighted.

==2 THE FINITE CALCULUS METHOD IN MECHANICS==

Let us consider a 1D space-time domain <math display="inline">\bar \Omega : \Omega \times t</math>, where <math display="inline">\Omega </math> is the domain length and <math display="inline">t</math> is the time.

We can now invoke any of the classical balance laws in mechanics, to be satisfied within a space-time slab of finite size <math display="inline">[x-d,x]\times [t-T ,t] \subset \bar \Omega </math> where <math display="inline">d</math> is the length of the space balance domain and and <math display="inline">T</math> is a time increment defining the size of the balance domain in the time axis.  Using the standard procedure of expanding in Taylor series the changes of the chosen governing variables and retaining the linear terms in <math display="inline">d</math> and <math display="inline">T</math> gives the FIC  balance equations as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r- \frac{h}{2}{\partial r \over \partial x}-\frac{\delta }{2}{\partial r \over \partial t}=0\qquad \hbox{in } \bar \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">r</math> denotes the governing differential equations of the infinitessimal theory. In Eq.(1) <math display="inline">h</math> and <math display="inline">\delta </math> are respectively characteristic length and characteristic time parameters satisfying

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}-d \le h\le d\\ 0\le \delta \le T\end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

For an advection-diffusion thermal problem, Eq.(1) is the heat balance equation with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r:=-\rho c \left(\displaystyle {\partial \phi  \over \partial t}+v \displaystyle {\partial \phi  \over \partial x}\right)+ \displaystyle {\partial  \over \partial x}\left(k \displaystyle {\partial \phi  \over \partial x}\right)+ Q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where the <math display="inline">\phi </math> is the temperature, <math display="inline">\rho </math> is the density, <math display="inline">v</math> is the velocity, <math display="inline">k</math> and <math display="inline">c</math> are the thermal diffusivity and the specific heat  parameters, respectively and <math display="inline">Q</math> is the heat source term.

Let us consider now a 1D solid mechanics problem. In this case Eq.(1) is the force equilibrium equation with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r:=-\rho {\partial ^2 u\over \partial t^2} + {\partial N \over \partial x}+b </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where <math display="inline">u</math> is the displacement, <math display="inline">N</math> is the axial force and <math display="inline">b</math> is the axial load.

The derivation of Eq.(1) for a simple steady state problem is given in the Appendix. Details of the derivation of the FIC equations for the transient case can be found in <span id='citeF-6'></span><span id='citeF-7'></span>[[#cite-6|[6,7]]].

Note that as the dimensions of the balance slab tend to zero, also the values of the characteristic parameters <math display="inline">h</math> and <math display="inline">\delta </math> tend to zero. In the limit case Eq.(1) tends to

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r=0\qquad \hbox{in } \bar \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

which is the standard form of the governing differential equations of the infinitessimal theory.

Eq.(1) can be generalized for a multidimensional problem as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_i - \frac{h_j}{2}{\partial r_i \over \partial x_j}-\frac{\delta }{2}{\partial r_i \over \partial t}=0 \quad  i=1,n_b\quad ;\quad j=1,n_d </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where <math display="inline">n_b</math> and <math display="inline">n_d</math> are the number of balance equations and the number of space dimensions, respectively.

A typical example of Eq.(6) are the FIC equations for an incompressible fluid written for the 3D case as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} - {h_j\over 2}{\partial r_{m_i} \over \partial x_j} - {\delta \over 2} {\partial r_{m_i} \over \partial t}=0</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|-
| style="text-align: center;" | <math> r_d -  {h_j\over 2}{\partial r_d \over \partial x_j}- {\delta \over 2} {\partial r_d \over \partial t}=0 \quad \quad i,j=1,2,3 </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}
|}

In above <math display="inline">r_{m_i}</math> and <math display="inline">r_d</math> are the ith momentum equation and the mass balance equation given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{m_i} := - \rho \left({\partial v_i \over \partial t}+v_j{\partial v_i \over \partial x_j}\right)+{\partial \sigma _{ij} \over \partial x_j}+b_i</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|-
| style="text-align: center;" | <math> r_d := {\partial v_i \over \partial x_i} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}
|}

where <math display="inline">v_i</math> are the velocities, <math display="inline">\sigma _{ij}=\tau _{ij}-p\delta _{ij}</math> are the total stresses, <math display="inline">p</math> is the pressure, <math display="inline">\tau _{ij}</math> are the viscous stresses and <math display="inline">b_i</math> are the body forces in the fluid <span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-3|[3,4]]].

Note that Eq.(8) is the result of applying the FIC concepts to the equations expressing the balance of mass over a space-time domain of finite size <span id='citeF-2'></span><span id='citeF-4'></span>[[#cite-2|[2,4]]].

For a general solid mechanics problem the governing FIC equations stating the equilibrium of forces will be identical to Eq.(7) where <math display="inline">r_{m_i}</math> is given by Eq.(9), dropping now the convective terms <span id='citeF-8'></span><span id='citeF-9'></span>[[#cite-8|[8,9]]].

===2.1 Boundary and initial conditions===

The Dirichlet and initial conditions are the same as in the standard infinitessimal theory. The Neumann boundary condition expressing balance of fluxes across the boundary <math display="inline">\Gamma _q</math> is however modified in the FIC formulation introducing, for consistency, the concept of balance over a domain of finite size adjacent to the Neumann boundary <math display="inline">\Gamma _q</math> <span id='citeF-2'></span>[[#cite-2|[2]]]. The new boundary condition is written generally as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_{n_i}-\bar q_{n_i} - {1\over 2} h_j n_j r_i =0\quad  \hbox{on }\Gamma _q \quad i=1,n_b\quad ;\quad j=1,n_d </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

where <math display="inline">q_{n_i}</math> is the normal “flow” across the boundary <math display="inline">\Gamma _q</math>, <math display="inline">\bar q_{n_i}</math> is the prescribed flow value and <math display="inline">n_i</math> are the components of the unit normal to the boundary.

For example, in a 3D advection-diffusion problem, the boundary condition for the “diffusive” flux is

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>n_j k_j {\partial \phi  \over \partial x_j}+\bar q_n - {1\over 2} h_j n_j r =0 \quad \quad j=1,2,3 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

For 3D fluid or solid mechanics problems, Eq.(11) would be

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma _{ij} n_j - \bar t_i - {1\over 2} h_j n_j r_{m_i} =0 \quad \quad i,j=1,2,3 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">\bar t_i</math> are the prescribed boundary tractions and <math display="inline">r_{m_i}</math> is given by Eq.(9).

Applications of the FIC method to problems in fluid and solid mechanics using the finite element method and the meshless finite point method can be found in <span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-7'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><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><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-2|[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]]].

==3 WHAT DO THE CHARACTERISTIC PARAMETERS MEAN?==

The characteristic parameters in the space and time dimension (<math display="inline">h_i</math> and <math display="inline">\delta </math>) can be interpreted as free ''intrinsic parameters'' giving the “exact” expression of the balance equations (up to first order terms) in a space-time domain of finite size.

Let us consider now the discretized solution of the modified governing equations. For simplicity we will focuss here in the simplest scalar 1D problem. The variable <math display="inline">\phi </math> is approximated as <math display="inline">\phi \simeq \hat \phi </math> where <math display="inline">\hat \phi </math> denotes the approximated solution. The values of <math display="inline">\hat \phi </math> are now expressed in terms of a finite set of parameters using any discretization procedure (finite elements, finite diferences, etc.). The discretized FIC governing equation (1) would read now

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat r- \frac{h}{2}{\partial \hat r \over \partial x}-\frac{\delta }{2}{\partial \hat r \over \partial t}=r_{\bar \Omega }\quad , \qquad \hbox{in } \bar \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">\hat r:=r (\hat \phi )</math> and <math display="inline">r_{\bar \Omega }</math> is the residual of the discretized FIC equations.

The meaning of the characteristic parameters <math display="inline">h</math> and <math display="inline">\delta </math> in the discretized equation (14) changes  (although the same simbols than in Eq.(1) have been kept). The <math display="inline">h</math> and <math display="inline">\delta </math> parameters become now ''of the order of magnitude of the discrete domain'' where the balance laws are satisfied. In practice, this means that

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}-l \le h\le l\\ 0\le \delta \le \Delta t\end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">l</math> is the grid dimension in the space discretization (i.e. the element size in the FEM or the cell size in the FDM) and <math display="inline">\Delta t</math> is the time step used to solve the transient problem.

The values of the characteristic parameters  can be found now in order to obtain an ''a correct numerical solution''. The meaning of “correct solution” must be obviously properly defined. Ideally, this would be a solution giving “exact” values at a discrete number of points (i.e. the nodes in a FE mesh). This is infortunatelly impossible for practical problems (with the exception of very simple 1D and 2D cases) and, in practice, <math display="inline">h</math> and <math display="inline">\delta </math> are computed making use of some chosen optimality rule, such as ensuring that the error of the numerical solution diminishes for appropriate values of the characteristic parameters <span id='citeF-2'></span><span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-2|[2,12,13]]]. This suffices in practice to obtain physically sound (stable) numerical results for any range of the physical parameters of the problem, always within the limitation of the discretization method chosen.

The concept of discretization of the FIC governing equations is explained in more detail with an example in the next section.

==4 FINITE ELEMENT DISCRETIZATION OF THE FIC EQUATIONS. AN EXAMPLE: THE ADVECTIVE-DIFFUSIVE PROBLEM==

Let us consider the FIC governing equations for the steady-state advective-diffusive problem defined in vector form for the multidimensional case as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r - {1\over 2}\mathbf{h}^T {\boldsymbol \nabla }r=0\quad \hbox{in } \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi -\bar \phi =0 \quad \hbox{on } \Gamma _\phi </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17a)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{n}^T \mathbf{D} {\boldsymbol \nabla }\phi + \bar{\mathbf{q}} -  {1\over 2}\mathbf{h}^T\mathbf{n} r  =0 \quad \hbox{on } \Gamma _q</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17b)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r:= - \mathbf{v}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla } ^T \mathbf{D} {\boldsymbol \nabla }\phi + Q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

In above equations <math display="inline">\boldsymbol h</math> is the characteristic length vector, <math display="inline">{\boldsymbol \nabla }</math> is the gradient vector, <math display="inline">\mathbf{D}</math> is the diffusivity matrix, <math display="inline">\mathbf{n}</math> is the normal vector and <math display="inline">\mathbf{v}</math> is the velocity vector. As usual <math display="inline">\Gamma _\phi </math> is the Dirichlet boundary where the unknown is prescribed to a fixed value <math display="inline">\bar \phi </math>.

A finite element interpolation of the unknown <math display="inline">\phi </math> can be written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi \simeq \hat \phi =\sum N_i \hat \phi _i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

where <math display="inline">N_i</math> are the shape functions and <math display="inline">\hat \phi _i</math> are the nodal values of the approximate function <math display="inline">\hat \phi </math> <span id='citeF-10'></span>[[#cite-10|[10]]].

Application of the Galerkin FE method to Equations (12)-(13b) gives, after integrating by parts the term <math display="inline">{\boldsymbol \nabla }r</math>  (and neglecting the space derivatives of <math display="inline">\mathbf{h}</math>) <span id='citeF-2'></span><span id='citeF-10'></span>[[#cite-2|[2,10]]]

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega N_i \hat r d\Omega - \int _{\Gamma _q} N_i (\mathbf{n}^T \mathbf{D} {\boldsymbol \nabla }\hat \phi +\bar q_n)d\Gamma + \sum \limits _e {1\over 2}  \int _{\Omega ^e} \mathbf{h}^T{\boldsymbol \nabla } N_i \hat r d\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

The last integral in Equation (20) has been expressed as the sum of the element contributions to allow for interelement discontinuities in the term <math display="inline">{\boldsymbol \nabla }\hat r</math>, where <math display="inline">\hat r = r (\hat \phi )</math> is the residual of the FE solution of the infinitesimal equations.

Note that the residual terms have disappeared from the Neumann boundary <math display="inline">\Gamma _q</math>. This is due to the consistency of the FIC terms in Equation (17b).

Integrating by parts the diffusive terms in the first integral of (20) leads to

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \left[N_i \mathbf{v}^T \nabla \hat \phi + {\boldsymbol \nabla }^T N_i \mathbf{D}{\boldsymbol \nabla }\hat \phi \right]d\Omega - \sum \limits _e {1\over 2} \int _{\Omega ^e} \mathbf{h}^T{\boldsymbol \nabla } N_i \hat r d\Omega - \int _\Omega N_i Q d\Omega +\int _{\Gamma _q} N_i \bar q_n d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

In matrix form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{K}\mathbf{a}=\mathbf{f} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

Matrix <math display="inline">\mathbf{K}</math> and vector <math display="inline">\mathbf{f}</math> are assembled from the element contributions given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K_{ij}^e = \int _{\Omega ^e} [N_i\mathbf{v}^T {\boldsymbol \nabla }N_j + {\boldsymbol \nabla }^T N_i (\mathbf{D} + {1\over 2}\mathbf{ h}^T\mathbf{v}){\boldsymbol \nabla } N_j]d\Omega - {1\over 2} \int _{\Omega ^e} \mathbf{h}^T {\boldsymbol \nabla } N_i {\boldsymbol \nabla } (\mathbf{D} {\boldsymbol \nabla } N_j)d\Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_i^e = \int _{\Omega ^e} [N_i + {1\over 2} \mathbf{h}^T{\boldsymbol \nabla } N_i] Qd\Omega - \int _{\Gamma _q^e}N_i \bar q_n d\Gamma  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

Note that the method introduces in <math display="inline">\boldsymbol K</math> an additional diffusivity matrix given by <math display="inline">{1\over 2} \mathbf{h}^T \mathbf{v}</math>. Also the second integral of Equation (23) vanishes for linear approximations. The same happens with the second term of the first integral of Equation (24) when <math display="inline">N_i</math> is linear and <math display="inline">Q</math> is constant. The evaluation of these integrals is mandatory in any other case.

===4.1 The discrete 1D advective-diffusive problem===

The role of the stabilization parameter can be better clarified with a simple 1D steady-state advection-diffusion problem governed by the FIC equation

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r - {h\over 2} {dr\over dx}=0 \quad \hbox{in } \quad 0\le x \le L </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

where <math display="inline">L</math> is the length of the analysis domain and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r = -v {d\phi \over dx}+k {d^2\phi \over dx^2}+Q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

The form of the element stiffness matrix and the element flux vector for the 1D case are obtained by simplification of Eqs.(23) and (24) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K_{ij}^e =\int _{l^e} \left[v N_i {\partial N_j \over \partial x} + \left(k+{vh\over 2}\right){\partial N_i \over \partial x} {\partial N_j \over \partial x}\right]dx - {1\over 2} \int _{l^e} h  {\partial N_i \over \partial x}k \frac{\partial ^2 N_j}{\partial x^2}dx</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27a)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_i^e = \int _{l^e} \left(N_i + {h\over 2}  {\partial N_i \over \partial x}\right)Q \,dx - \bar q_L</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27b)
|}

For simplicity we will assume <math display="inline">Q=0</math> <math display="inline">\phi=0</math> at <math display="inline">x=0</math> and <math display="inline">\phi = \bar \phi </math> at <math display="inline">x=L</math>.

We will solve the problem with the simplest  mesh of two linear finite element of dimensions <math display="inline">l^e=L/2</math>. The stencil for this problem is

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{v\over 2} (\phi _3 -\phi _1) + \left(k+{vh\over 2}\right)\left({-\phi _3 + 2\phi _2-\phi _1)\over l^e}\right)=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

The solution for the central node is obtained making <math display="inline">\phi _1=0</math> and <math display="inline">\phi _3=\bar \phi </math> in Eq.(28) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi _2 =\left[2(1+\alpha \gamma )\right]^{-1} (1-\gamma + \alpha \gamma )  \bar \phi  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

where <math display="inline">\gamma =\displaystyle{vl^e\over 2k}</math> is the element Peclet number and <math display="inline">\alpha =\displaystyle{h\over l^e}</math> is a dimensionless characteristic parameter.

We note that for <math display="inline">\alpha =0</math> (i.e. the discrete solution of the standard infinitessimal equations) <math display="inline">\phi _2 = {1\over 2} (1-\gamma ) \bar \phi </math> which gives non physical (unstable) negative values of <math display="inline">\phi _2</math> for <math display="inline">\gamma > 1</math>. A critical value of <math display="inline">\alpha </math> can be deduced from Eq.(29) ensuring that <math display="inline">\phi _2 \ge 0</math> for all range of the parameters <math display="inline">u,k</math> and <math display="inline">l</math>. This gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha \ge 1 - {1\over \gamma } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

which is the well known expression of the ''critical stabilization parameter'' for this problem <span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-10|[10,11]]].

The value of <math display="inline">\alpha </math> given by Eq.(30) allow us to obtain the best possible physically meaningful solution for the problem. Indeed for this simple case we could aim to obtain the ''exact analytical solution'' at node <math display="inline">\phi _2</math>. This is simply granted is <math display="inline">\alpha </math> is taken as <math display="inline">\alpha =\coth \gamma -  1/\gamma </math> which can be again deduced from Eq.(28) <span id='citeF-6'></span><span id='citeF-10'></span>[[#cite-6|[6,10]]]. Note that the values of this expression for <math display="inline">\alpha </math> and those of Eq.(30) are quite similar for <math display="inline">\gamma > 2</math> <span id='citeF-10'></span>[[#cite-10|[10]]].

The merit of the expression of <math display="inline">\alpha </math> of Eq.(30) is that it has been deduced ''with no previous knowledge of the exact solution''. Indeed, the approximated solution within the element adjacent to the boundary at <math display="inline">x=L</math> will differ substantially from the exact one: the numerical solution is linear whereas the analytical solution gives a ''boundary layer'' adjacent to <math display="inline">x=L</math>. The stabilized FIC numerical solution can however be taken now as the starting point to obtain a refined solution giving a more accurate description of the physical problem within the element where the boundary layer is located. This can be obtained by mesh refinement or either by enriching the approximation within that element in an adaptive manner. Here a number of possibilities can be explored using mesh refinement, hierarchical approximations, <math display="inline">p</math>-refinement or enrichment techniques, etc. Some of these options are discussed in a next section.

===4.2 The characteristic length vector===

From the FIC expressions in Section 4 we note that in multidimensional problems the space stabilization parameters can be grouped in a vector <math display="inline"> \boldsymbol h</math>, where <math display="inline"> \boldsymbol h</math> has the space dimensions of the problem (i.e. <math display="inline">\mathbf{h} =[h_1,h_2,h_3]^T</math> for 3D situations). Computation of vector <math display="inline"> \boldsymbol h</math> is a crucial step as the quality of the stabilized solution depends on the choice of the module and direction of <math display="inline"> \boldsymbol h</math>.

We could, for instance,  assume that the direction of vector <math display="inline">\boldsymbol h</math> is ''parallel'' to the velocity <math display="inline">\boldsymbol v</math>, i.e. <math display="inline">\mathbf{h}= h{\mathbf{v}\over \vert \mathbf{v}\vert }</math> where <math display="inline">h</math> is a characteristic length. Under these conditions, Equation (20) reads

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega N_i \hat r d\Omega - \int _{\Gamma _q} N_i (\mathbf{n}^T \mathbf{D} {\boldsymbol \nabla }\hat \phi +\bar q_n)d\Omega{+} \sum \limits _e \int _{\Omega ^e} {h\over 2\vert \mathbf{v}\vert  }   \mathbf{v}^T{\boldsymbol \nabla } N_i \hat rd\Omega =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

Equation (31) coincides precisely with the so called Streamline-Upwind-Petrov-Galerkin (SUPG) method <span id='citeF-2'></span><span id='citeF-10'></span><span id='citeF-20'></span>[[#cite-2|[2,10,20]]]. The ratio <math display="inline">\displaystyle{h\over 2\vert \mathbf{v}\vert  }</math> has dimensions of time and it is termed element ''intrinsic time'' parameter <math display="inline">\tau </math>.

It is important to note that the SUPG expression is a ''particular case'' of the more general FIC formulation. This explains the limitations of the SUPG method to provide stabilized numerical results in the vicinity of sharp gradients of the solution transverse to the flow direction. In general, the adequate direction of <math display="inline">\mathbf{ h}</math> is not coincident with that of <math display="inline">\mathbf{v}</math> and the components of <math display="inline">\mathbf{h}</math> introduce the necessary stabilization along the streamlines and the transverse directions to the flow. In this manner, the FIC method reproduces the best-features of the so-called stabilized discontinuity-capturing schemes <span id='citeF-5'></span><span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-5|[5,12,13]]].

The computation of the characteristic vector <math display="inline"> \boldsymbol h</math> is a topioc which fails outside the scope of this paper. We will just mention that in practice it is convenient to split <math display="inline"> \boldsymbol h</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{h} = \mathbf{h}_s + \mathbf{h}_t </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

where <math display="inline">\mathbf{h}_s</math> and <math display="inline">\mathbf{h}_t</math> are called streamline and transverse characteristic lenght vectors. A useful option is to take

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{h}_s= h_s {\mathbf{v}\over \vert \mathbf{v}\vert } \qquad \hbox{and}\quad \mathbf{ h}_t=h_t{{\boldsymbol \nabla }\phi \over \vert {\boldsymbol \nabla }\phi \vert } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

The characteristic lengths <math display="inline">{h}_s</math> and <math display="inline">{h}_t</math> are computed at element level by imposing that the residual of the numerical solution diminishes over each element <span id='citeF-2'></span><span id='citeF-5'></span><span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-2|[2,5,12,13]]]. In practice, these lengths can be computed for each element as the maximum value of the scalar product between the element sides and the velocity vector (for <math display="inline">h_s</math>) and the gradient vector (for <math display="inline">h_t</math>) <span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span>[[#cite-14|[14,15,16,17]]].

Note that the gradient of <math display="inline">\phi </math> in the expression of <math display="inline">\mathbf{h}_t</math> introduces a non linearity in the computation process which is typical of discontinuity capturing techniques. A procedure to linealize the expression of the <math display="inline">\mathbf{h}_t</math> in element next to boundary layers by approximating <math display="inline">{\boldsymbol \nabla }\phi </math> by the normal vector to the boundary is presented in <span id='citeF-5'></span>[[#cite-5|[5]]].

==5 INTERPRETATION OF THE DISCRETE SOLUTION OF THE FIC EQUATIONS==

Let us consider the solution of a physical problem in a space domain <math display="inline">\Omega </math>, governed by a differential equation <math display="inline">r(\phi ) =0</math> in <math display="inline">\Omega </math> with the corresponding boundary conditions. The “exact” (analytical) solution of the problem will be a function giving the sought distribution of <math display="inline">\phi </math> for any value of the geometrical and physical parameters of the problem. Obviously, since the analytical solution is difficult to find (practically impossible for real situations), an approximate numerical solution is found <math display="inline">\phi \simeq \hat \phi </math> by solving the problem <math display="inline">\hat r =0</math>, with <math display="inline">\hat r = r(\hat \phi )</math>, using a particular discretization method (such as the FEM). The distribution of <math display="inline">\phi </math> in <math display="inline">\Omega </math> is now obtained for specific values of the geometrical and physical parameters. The accuracy of the numerical solution  depends on the discretization parameters, such as the number of elements and the approximating functions chosen in the FEM. Figure 1 shows a schematic representation of the distribution of <math display="inline">\hat \phi </math> along a line for different discretizations <math display="inline">M_1,M_2,\cdots , M_n</math> where <math display="inline">M_1</math> and <math display="inline">M_n</math> are the coarser and finer meshes, respectively. Obviously for <math display="inline">n</math> sufficiently large a good approximation of <math display="inline">\phi </math> will be obtained and for <math display="inline">M_\infty </math> the numerical solution <math display="inline">\hat \phi </math> will coincide with the “exact” (and probably unreachable) analytical solution <math display="inline">\phi </math> at all points. Indeed in some problems the <math display="inline">M_\infty </math> solution can be found by a “clever” choice of the discretization parameters.

An unstable solution will occur when for some (typically coarse) discretizations, the numerical solution provides non physical or very unaccurate values of <math display="inline"> \hat \phi </math>. A situation of this kind is represented by curves <math display="inline">M_1</math> and <math display="inline">M_2</math> of the left hand side of Figure 1. These unstabilities would disappear by an appropriate mesh refinement  (curves <math display="inline">M_3,M_4 \cdots </math> in Figure 1) at the obvious increase of the computational cost.

In the FIC formulation the starting point are the modified differential equations of the problem and the corresponding modified boundary conditions as previously described. The modified equations are not useful to find an analytical solution, <math display="inline">\phi (x)</math>, for the physical problem. However, the numerical solution of the FIC equation can be readily found. Moreover, by adequately choosing the values of the characteristic length parameter <math display="inline">h</math>, the numerical solution of the FIC equations will be always stable (physically sound) for any discretization level chosen.

This process is schematically represented in Figure 1 where it is shown that the numerical oscillations  for the coarser discretizations <math display="inline">M_1</math> and <math display="inline">M_2</math> dissapear when using the FIC procedure.

We can therefore conclude the FIC approach allows us to obtain a ''better numerical solution for a given discretization''. Indeed, as in the standard infinitesimal case, the choice of <math display="inline">M_\infty </math> will yield the (unreachable) exact analytical solution and this ensures the consistency  of the method.

<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:Draft_Samper_937954379-Fig13con215.png|456px|Schematic representation of the numerical solution of a physical problem using standard infinitesimal calculus and finite calculus.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Schematic representation of the numerical solution of a physical problem using standard infinitesimal calculus and finite calculus.
|}

==6 MULTISCALE ADAPTIVE FINITE ELEMENT ANALYSIS USING THE FIC FORMULATION==

===6.1 Adaptive solution by updating the stability parameters===

The FIC formulation can be applied in conjunction with any numerical method (say the finite element method) in order to obtain a “correct” solution within the limit of the discretization chosen. The accuracy of the finite element solution depends greatly of the values chosen for the stabilization parameters in the FIC formulation. Indeed if these values are taken equal to zero, an unstable non physical solution is obtained in many cases. The space and time stabilization parameters <math display="inline">h_i</math> and <math display="inline">\delta </math> in the FIC formulation can be given adequate values such that the desired accuracy in the numerical solution is achieved. In <span id='citeF-2'></span><span id='citeF-5'></span><span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-2|[2,5,12,13]]] a procedure is presented to adaptively compute the values of the stabilization parameters so that the element error, measured in terms of the average residual of the governing equations over the element, diminishes until it reaches a prescribed value. Other procedures to compute the <math display="inline">h_i</math> and <math display="inline">\delta </math> parameters can be devised, all aiming to obtain the “best” possible nodal solution for a given mesh.

Once this solution has been obtained, the problem of accurately reproducing sharp gradients in selected zones remains. Some methods aiming to solve this problem are proposed next.

The nodal finite element values computed with the stabilized FIC formulation can be taken as the starting point for obtaining an enhanced solution in zones where high gradients of the sought variables occur. Three options will be considered here: mesh adaptivity, hierarchical shape functions and enriched interpolations incorporating the shape of the analytical solution.

===6.2 Mesh refinement===

This is the more standard procedure, although it is still unpractical for many practical 3D problems involving complex geometries. The mesh adaption process can be based on mesh regeneration or mesh enrichment procedures. The adaptivity strategy can use global or local error estimation based on adequate error norms. Details of  mesh adaptivity procedures based on different error measures can be found elsewhere and will not be discussed here <span id='citeF-10'></span><span id='citeF-21'></span><span id='citeF-22'></span>[[#cite-10|[10,21,22]]]. Given the intrinsic difficulties for modifying the meshes for 3D problems we will favour the two procedures described next versus the mesh adaptivity process.

===6.3 Hierarchical interpolation===

The nodal values provided by the FIC method can be considered “correct” (in the sense that they are physically sound) whithin the limitations of the mesh discretization chosen. A procedure to obtain an enhanced approximation within a particular element, or a patch of elements, is to introduce hierarchical interpolations in the selected elements. The hierarchical interpolation can be simply written for an individual element

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi = \sum \limits _{i=1}^n N_i \phi _i + \sum \limits _{j=n+1}^m N_j^h a_j^h </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

where <math display="inline">N_i</math> denote the standard shape functions, <math display="inline">\phi _i</math> are the nodal values, <math display="inline">N_j^h</math> are the hierarchical shape functions and <math display="inline">a_j^h</math> are the hierarchical values. As usual, the functions <math display="inline">N_j^h</math> should vanish along the element edges. The simplest hierarchical interpolations (satisfying orthogonality conditions) can be based on Legendre polynomials or on standard trigonometric functions <span id='citeF-10'></span>[[#cite-10|[10]]]. More advanced hierarchical finite element interpolations based on the concept of the partition of unity can be effectively used <span id='citeF-23'></span>[[#cite-23|[23]]].

A good approximation of the hierarchical variables can be obtained taking the nodal values obtained with the FIC numerical solution as ''fixed''. A simple set of equations needs to be solved at each element level only for the hierarchical variables <math display="inline"> a_j^h</math>. Indeed, this computation can be performed within a more general iterative procedure, where the nodal variables and the hierarchical values are updated in an staggered manner.

===6.4 Enrichment interpolation using an approximation to the analytical shape of the solution===

An effective procedure to obtain an enhanced solution starting with the FIC nodal values is to use a local interpolation over each element injecting into the interpolating functions a good approximation to the shape of the sought solution. This can be effectively carried out by using as the new approximation the free-space solution of the homogeneous and constant-coefficient version of the standard governing partial differential equations. A procedure of this kind has been used by Farhat ''et al.'' <span id='citeF-24'></span>[[#cite-24|[24]]] who added the enriched field to the initial polynomial field. The resulting field is not continuous in this case and continuity is enforced by means of Lagrange multipliers.

In this work it is proposed to apply the “analytical” enrichment in order to obtain an enhanced solution within the selected elements where high gradients occur, while keeping the nodal values given by the FIC solution fixed.

Indeed the three adaptive procedures above suggested can be combined and the efficiency of the different alternatives should be investigated and tested in practice.

The concepts explained are clarified next with two simple examples of application.

==7 EXAMPLES==

===7.1 1D advective-diffusive problem===

We will consider the simplest 1D advective-diffusive problem governed by Eqs.(25) and (26) with boundary conditions <math display="inline">\phi =0</math> at <math display="inline">x=0</math>, <math display="inline">\phi =\bar \phi </math> at <math display="inline">x=L</math>. Again for simplicity we will consider a mesh of two linear 1D elements.

Figure 2 shows the unstable numerical solution obtained for the case of <math display="inline">\alpha =0</math> (Galerkin) and the “correct” FIC solution for a critical value of <math display="inline">\alpha =(1- {1\over \gamma }) (1+\varepsilon )</math> with <math display="inline">\varepsilon =10^{-10}</math>. Note that the numerical solution, although nodally correct, is far from the “exact”  analytical distribution giving a boundary layer next to the end node at <math display="inline">x=L</math>.

<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:Draft_Samper_937954379-Fig2.png|600px|1D advection-diffusion problem solved with two linear elements (γ=50). Unstable solution (- - -) obtained for α=0 (Galerkin). Stabilized solution (&#8211;-) obtained with the FIC method for α= (1+ 10⁻¹⁰)(1- 1\over γ)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' 1D advection-diffusion problem solved with two linear elements (<math>\gamma =50</math>). Unstable solution (- - -) obtained for <math>\alpha =0</math> (Galerkin). Stabilized solution (&#8211;-) obtained with the FIC method for <math>\alpha = (1+ 10^{-10})(1- {1\over \gamma })</math>
|}

We will consider next a hierarchical approximation of the solution over the element 2&#8211;3 using Eq.(34) with hierarchical interpolation functions defined as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>N_j^h =\sin j {\pi \over 2} (1+\xi ) \quad \hbox{with } -1 \le \xi \le 1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

Note that <math display="inline">N_j^h</math> vanishes for the element ends located at <math display="inline">\xi =\pm 1</math>. Figure 3 shows the approximation of the numerical solution over the element 2&#8211;3 obtained for <math display="inline">m=3</math> ''keeping the nodal values given by the original FIC solution fixed''. The improvement in the solution is noticeable. An almost identical solution is obtained if the value for <math display="inline">\phi _2</math> is recomputed until convergence is achieved. A similar improvement is obtained using Legendre polynomials up to a quartic order.

A considerable better solution over the element 2&#8211;3 can be obtained by using for approximating the unknown within this element the following interpolation

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi = \sum \limits _{i=1}^2 P_i (\xi )  \phi _i^{(e)} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

with <math display="inline">P_1 = \displaystyle{e^{\gamma \xi }-e^{\gamma }\over e^{-\gamma } - e^{\gamma }}\quad ,\quad P_2=\displaystyle{e^{-\gamma } - e^{\gamma \xi }\over e^{-\gamma } - e^{\gamma }}</math>.

The new shape functions <math display="inline">P_1</math> and <math display="inline">P_2</math> are obtained from the analytical solution of the homogeneous differential equation (<math display="inline">r=0</math>) for constant values of <math display="inline">u</math> and <math display="inline">k</math>. Note that <math display="inline">P_i (\xi _j ) =1</math> for <math display="inline">i=j</math> and <math display="inline">=0</math> for <math display="inline">i\not = j</math>.

Application of the approximation (36) over the element 2&#8211;3 with <math display="inline">\phi _2^{(e)} =\bar \phi </math> and <math display="inline">\phi _1^{(e)} = \phi _2</math>, where <math display="inline">\phi _2</math> is the nodal value at node 2 obtained with the FIC formulation, reproduces accurately the boundary layer within this element, as expected (see Figure 3). Recall that Eq.(36) has been applied “a posteriori”, once the nodal value at node 2 has been found using the original linear interpolation and the FIC method with the critical value of <math display="inline">\alpha </math> above given.

<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:Draft_Samper_937954379-Fig3.png|600px|1D advection-diffusion problem solved with two linear elements. Initial stabilized FIC solution (curve A). Enhanced solutions over element 2&#8211;3 using a hierarchical trigonometric interpolation (curve B) and an exponential interpolation (curve C)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' 1D advection-diffusion problem solved with two linear elements. Initial stabilized FIC solution (curve A). Enhanced solutions over element 2&#8211;3 using a hierarchical trigonometric interpolation (curve B) and an exponential interpolation (curve C)
|}

This simple example shows that the nodal values provided by the FIC method are a good starting point for obtaining an enhanced numerical solution in the vecinity of sharp gradients. The results from this simple example indicate that it suffices with a straightforward “a posteriori” analysis to obtain a better solution in these regions using either a hierarchical interpolation or (ideally) an approximation enriched with the shape of the solution of the homogeneous equation for constant coefficients.

===7.2 2D stationary convection-diffusion problem with diagonal velocity, zero source and uniform Dirichlet conditions===

<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:Draft_Samper_937954379-Figura1a.png|312px|]]
|[[Image:Draft_Samper_937954379-Figura1b.png|300px|]]
|-
|[[Image:Draft_Samper_937954379-Figura1c.png|390px|]]
|[[Image:Draft_Samper_937954379-Figura1d.png|438px|Convection-diffusion problem with uniform Dirichlet conditions. One step solution. Distribution of ϕ along the lines AA' and BB']]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' Convection-diffusion problem with uniform Dirichlet conditions. One step solution. Distribution of <math>\phi </math> along the lines <math>AA'</math> and <math>BB'</math>
|}

The steady-state convection-diffusion equation is solved in a domain of unit size  (Figure 4) with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k_x = k_y = 1\quad , \quad \mathbf{v}=10^{10} [1,1]^T\quad ,\quad Q=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

The following Dirichlet conditions are assumed

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi =0 \quad \hbox{on the lines }\quad x=0 \quad \hbox{and}\quad y=0</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi =100 \quad \hbox{on the lines}\quad x=1 \quad \hbox{and}\quad y=1</math>
|}
|}

The expected solution is a uniform distribution of <math display="inline">\phi =0</math> over the domain except in the vecinity of the boundaries <math display="inline">x=1</math> and <math display="inline">y=1</math> where a boundary layer is formed.

The domain is discretized with a uniform mesh of 400 four node quadrilaterals (Figure 4).

The solution has been obtained in one single step using the value of <math display="inline">\mathbf{h}</math> given by Eqs.(32) and (33) approximating the gradient of <math display="inline">\phi </math> in the elements next to the outflow boundaries by the normal vector. For details see <span id='citeF-5'></span>[[#cite-5|[5]]]. Figure 4 shows that the sharp gradients expected are perfectly captured over the boundary elements without any oscillation.

An enhanced solution over the elements adjacent to the outflow boundaries can now be obtained using hierarchical interpolations or exponential interpolations similarly as described for the 1D example in the previous section. The efficiency of these procedures is similar to that shown in Figure 4 for the 1D case giving in both cases an enhanced reproduction of the boundary layer within the elements adjacent to the lines <math display="inline">x=1</math> and <math display="inline">y=1</math>.

==8 POSSIBILITIES OF FINITE CALCULUS FOR STRAIN LOCALIZATION==

The FIC method introduces naturally higher order derivative terms of the displacements in the equilibrium equations. These terms resemble those introduced by the Cosserat model <span id='citeF-25'></span>[[#cite-25|[25]]] and the so called “non local” constitutive models <span id='citeF-25'></span><span id='citeF-26'></span><span id='citeF-27'></span><span id='citeF-28'></span><span id='citeF-29'></span>[[#cite-25|[25,26,27,28,29]]]. These models are typically used in order to preserve the ellipticity of the solid mechanics equations in the presence of localized high displacement gradient zones, such as shear bands and fracture lines.

For instance, the FIC governing equations for a simple 1D bar can be written (in absence of external body forces) as <span id='citeF-5'></span>[[#cite-5|[5]]]

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\partial N \over \partial x}-{h\over 2}{\partial ^2N \over \partial x^2}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

where <math display="inline">N</math> is the axial force.

The relationship between the axial force and the elongation <math display="inline">\varepsilon =\displaystyle{du\over dx}</math> in the softening branch can be written in incremental form as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta N =- \vert E_T\vert \delta \varepsilon = - \vert E_T\vert  {d(\delta u)\over dx} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

where <math display="inline">E_T</math> is the tangent modulus of the material and <math display="inline">\delta u</math> is the displacement increment.

Substituting Equation (39) into (38) gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{d^2 (\delta u)\over dx^2} - {h\over 2} {d^3 (\delta u)\over dx^3}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

Equation (40) resembles the governing equations derived by Lasry and Belytschko <span id='citeF-26'></span>[[#cite-26|[26]]] using localization limiters in the stress-strain relationship and by Schreier and Chen <span id='citeF-27'></span>[[#cite-27|[27]]] using gradient-dependent plasticity constitutive models.

Many similarities can be found between the governing equations derived by the FIC method and those resulting from the non local constitutive models <span id='citeF-30'></span>[[#cite-30|[30]]]. Although in some situations the resulting governing differential equations are the same, the basic difference between the two approaches is that non local constitutive models aim to enhance the constitutive equation by adding new terms (typically strain gradient terms), whereas the FIC method defines a new equilibrium equation applicable to all the discrete scales while preserving the features of the constitutive equation. This opens a world of possibilities for the FIC equations in solid mechanics to solve strain localization problems using standard constitutive models.

==9 CONCLUSIONS==

The paper has presented some basic ideas to combine the finite calculus (FIC) method with adaptive finite element strategies in order to solve problems where sharp gradients of the solution occur in localized zones. The key argument of the procedures proposed is that the discretized FE solution of the FIC equations leads to nodal values which are “as good as possible” by virtue of the adequate selection of the stabilization parameters. The FIC solution can be used as the starting point to compute an “a posteriori” enhanced solution over selected elements using hierarchical or enriched approximations. Indeed, mesh adativity can be included here as a third alternative and the combination of the different procedures opens new areas of research.

The same ideas can be used to solve a wide class of multiscale problems in mechanics where the solution changes in several orders of magnitude in specific regions.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' Nottale, L. Relativity in all its states. (in French), Hachette, 1998.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Oñate E. Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. ''Comput. Methods Appl. Mech. Engrg.'' 1998; '''151:1-2''':233&#8211;267.

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  Oñate E,  García J. A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation. in ''Comput. Methods Appl. Mech. Engrg.'' 2001; '''191:6-7''':635-660.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  Oñate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. ''Comput. Methods Appl. Mech. Engrg.'' 2000; '''182:1&#8211;2''':355&#8211;370.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  Oñate E. Possibilities of finite calculus in computational mechanics. To be published in ''Int. J. Num. Meth. Engng.'', 2003.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Oñate E, Manzan M. Stabilization techniques for finite element analysis of convection diffusion problems. In ''Comput. Anal. of Heat Transfer'',  WIT Press, G. Comini and B. Sunden (Eds.), 2000.

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Oñate E, Manzan M. A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. ''Int. J. Num. Meth. Fluids'' 1999; '''31''':203&#8211;221.

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' Oñate E, Rojek J, Taylor RL, Zienkiewicz OC. Finite calculus formulation for analysis of incompressible solids using linear triangles and tetrahedra. Submitted to ''Int. J. Num. Meth. Engng.'', May 2002.

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  Oñate E, Rojek J, Taylor RL, Zienkiewicz OC. Linear triangles and tetrahedra for incompressible problems using a finite calculus formulation. In ''European Conference on Computational Mechanics (ECCM2001)'', Cracow, Poland, 26&#8211;29 June, 2001.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]''' Zienkiewicz OC, Taylor RL. ''The finite element method''. 5th Edition, 3 Volumes, Butterworth&#8211;Heinemann, 2000.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  Hirsch C. ''Numerical computation of internal and external flow'', J. Wiley, Vol. 1 1988, Vol. 2, 1990.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]''' Oñate E, García J, Idelsohn SR. Computation of the stabilization  parameter for the finite element solution of advective-diffusive problems. ''Int. J. Num. Meth. Fluids'' 1997;'''25''':1385&#8211;1407.

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]''' Oñate E, García J, Idelsohn SR.  An alpha-adaptive approach for stabilized finite element solution of advective-diffusive problems with sharp gradients. ''New Advances in Adaptive Comput. Methods in Mech.'', Elsevier, P. Ladeveze, J.T. Oden (Eds.), 1998.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  García J, Oñate E. An unstructured finite element solver for ship hydrodynamics. ''J. Appl. Mechanics'' 2002.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  Oñate E, García J, Bugeda G, Idelsohn SR. A general stabilized formulation for incompressible fluid flow using finite calculus and the FEM. In ''Towards a New Fluid Dynamics with its Challenges in Aeronautics'', J. Periaux, D. Champion, O. Pironneau and Ph. Thomas (Eds.), CIMNE, Barcelona, 2002.

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''   Oñate E, García J, Idelsohn SR. Ship Hydrodynamics. ''Encyclopedia of Computational Mechanics'', T. Hughes, R. de Borst and E. Stein (Eds.), J. Wiley, 2004 (to be published).

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''  Oñate E, Taylor RL, Zienkiewicz OC,   Rojek J. A residual correction method based on finite calculus. ''FEMClass of 42 Meeting'', 30&#8211;31 May, Ibiza, 2002, www.congress.cimne.upc.es/femclass42. To be published in Engineering Computation.

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  Oñate E, Idelsohn SR, Zienkiewicz OC, Taylor RL. A finite point method in computational mechanics. Applications to convective transport and fluid flow. ''Int. J. Num. Meth. Engng.'' 1996;'''39''':3839&#8211;3866.

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  Oñate E, Sacco C, Idelsohn S. A finite point method for incompressible flow problems. ''Computing and Visualization in Science'' 2000; '''2''':67&#8211;75.

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]''' Brooks A, Hughes TJR. Streamline upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. ''Comput. Methods Appl. Mech. Engrg.'' 1982; '''32''':199&#8211;259.

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  Ladeveze P., Pelle J. ''La maïtrise du calcul en mecanique linearise et non lineaire''. Hermes, 2001.

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  Oñate E, Bugeda G. A study of mesh adaptivity criteria in adaptive finite element computations. ''Engng. Computations'' 1993; '''10''':307&#8211;321.

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''  Taylor RL, Zienkiewicz OC, Oñate E. A hierarchical finite element method based on the partition of unity. ''Comput. Methods Appl. Mech. Engrg.'' 1998; '''152''':73&#8211;84

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  Farhat C, Harari I, Franca L. The discontinuous enrichment method. ''Comput. Methods Appl. Mech. Engrg.'' 2001; '''190''':6455&#8211;6479.

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  de Borst R. Simulation of strain localization: A reappraisal of the Cosserat continuum. ''Continuum Computation'' 1991; '''8''':317&#8211;332.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  Larsry D, Belytschko T. Localization limiters in transient problems. ''Int. J. Solids and Structures'' 1988; '''24''':581&#8211;597.

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  Schreier HL, Chen Z. One dimensional softening with localization. ''ASME J. Appl. Mech.'' 1980; '''53''':791&#8211;797.

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  de Borst R. Gradient-dependent plasticity: Formulation and algorithmic aspects. ''Int. J. Num. Meth. Engng'' 1992; '''35''':521&#8211;539.

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''  Peric D, Yu J, Owen DRJ. On error estimates and adaptivity in elastoplastic solids: Application to the numerical simulation of strain localization in classical and Cosserat continua. ''Int. J. Num. Meth. Engng.'' 1994; '''37''':1351&#8211;1379.

<div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''  Borja R. Possibilities of the FIC method for strain localization problems. Private communication, 2001.

==APPENDIX==

We will consider a convection-diffusion problem in a 1D domain <math display="inline">\Omega </math> of length <math display="inline">L</math>. The equation of balance of fluxes in a subdomain of size <math display="inline">d</math> belonging to <math display="inline">\Omega </math> (Figure A1) is written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_A - q_B=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

where <math display="inline">q_A</math> and <math display="inline">q_B</math> are the incoming and outgoing fluxes at points <math display="inline">A</math> and <math display="inline">B</math>, respectively. The flux <math display="inline">q</math> includes both convective and diffusive terms; i.e. <math display="inline">q=v\phi - k{d\phi \over dx}</math>, where <math display="inline">\phi </math> is the transported variable (i.e. the temperature in a thermal problem), <math display="inline">v</math> is the velocity and <math display="inline">k</math> is the diffusitivity of the material.

<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:Draft_Samper_937954379-majesus2.png|240px|Equilibrium of fluxes in a  space balance domain of finite size]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Equilibrium of fluxes in a  space balance domain of finite size
|}

Let us express now the fluxes <math display="inline">q_A</math> and <math display="inline">q_B</math> in terms of the flux at an arbitrary point <math display="inline">C</math> within the balance domain (Figure A.1). Expanding <math display="inline">q_A</math> and <math display="inline">q_B</math> in Taylor series around point <math display="inline">C</math> up to second order terms gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_A= q_C - d_1 \frac{dq}{d x}\vert _C+ \frac{d^2_1}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_1)\quad ,\quad  q_B= q_C + d_2 \frac{dq}{d x}\vert _C+\frac{d^2_2}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_2) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

Substituting Equation (A.2) into Equation (A.1) gives after simplification

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{dq}{dx}-\underline{\frac{h}{2} \frac{d^2q}{dx^2}}\frac{h}{2} \frac{d^2q}{dx^2}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

where <math display="inline">h=d_1-d_2</math> and all the derivatives are computed at point <math display="inline">C</math>.

Standard calculus theory assumes  that the domain <math display="inline">d</math> is of infinitesimal size and the resulting balance equation is simply <math display="inline">{dq\over dx}=0</math>. We will relax this assumption and allow the space balance domain to have a ''finite size''. The new balance equation (A.3) incorporates now the underlined term which introduces the ''characteristic length'' <math display="inline">h</math>. Obviously, accounting for higher order terms in Equation (A2) would lead to new terms in Equation (A.3) involving higher powers of <math display="inline">h</math>.

Distance <math display="inline">h</math> in Equation (A.3) can be interpreted as a free parameter depending on the location of point <math display="inline">C</math> within the balance domain. Note that <math display="inline">-d\le h\le d</math> and, hence, <math display="inline">h</math> can take a negative value. At the discrete solution level the domain <math display="inline">d</math> should be replaced by the balance domain around a node. This gives for an equal size discretization <math display="inline">-l^e\le h \le l^e</math> where <math display="inline">l^e</math> is the element or cell dimension. The fact that Equation (A.3) is the ''exact balance equation'' (up to second order terms) for any 1D domain of finite size and that the position of point <math display="inline">C</math> is arbitrary, can be used to derive numerical schemes with enhanced properties simply by computing the characteristic length parameter from an adequate “optimality” rule leading to an smaller error in the numerical solution.

Consider, for instance, Equation (A.3) applied to the 1D convection-diffusion problem. Neglecting third order derivatives of <math display="inline">\phi </math>, Equation (A.3) can be rewritten  in terms of <math display="inline">\phi </math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-v \frac{d\phi }{dx}+\left(k+\frac{v h}{2}\right)\frac{d^2\phi }{dx^2}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

We see clearly that the FIC method introduces ''naturally'' an additional diffusion term in the standard convection-diffusion equation. This is the basis of the popular “artificial diffusion” procedure <span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-10|[10,11]]] where the characteristic length <math display="inline">h</math> is typically expressed as a function of the cell or element dimension.

Equation (A.3) can be extended to account for source terms. The modified governing equation can then be  written in compact form as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r- \underline{\frac{h}{2} \frac{dr}{dx}}\frac{h}{2} \frac{dr}{dx}=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r : = -v \frac{d\phi }{dx}+ \frac{d}{dx}\left(k \frac{d\phi }{dx}\right)+ Q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

where <math display="inline">Q</math> is the external source. For consistency a “finite” form of the Neumann boundary condition should be used. This can be readily obtained by invoking balance of fluxes in a domain of finite size next to the boundary <math display="inline">\Gamma _q</math> where the external flux is prescribed to a value <math display="inline">q_p</math>. The modified Neumann boundary condition is <span id='citeF-2'></span>[[#cite-2|[2]]]

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k \frac{d\phi }{dx}+q_p - \underline{\frac{h}{2}r}\frac{h}{2}r=0 \quad \hbox{at } \Gamma _q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

The governing equations of the problem are completed with the standard Dirichlet condition prescribing the value of <math display="inline">\phi </math> at the boundary <math display="inline">\Gamma _\phi </math>.

The underlined terms in Equations (A.5) and (A.7) introduce the necessary stabilization in the discrete solution  using whatever numerical scheme.

The time dimension can be simply accounted for the FIC method by considering the balance equation in a space-time slab domain as described in <span id='citeF-2'></span><span id='citeF-6'></span><span id='citeF-7'></span>[[#cite-2|[2,6,7]]].