# An element-based formulation for ES-FEM and FS-FEM models for implementation in standard solid mechanics finite element codes for 2D and 3D static analysis

Daniele Colombo, Slah Drira, Ralf Frotscher, Manfred Staat

► To cite this version:

Daniele Colombo, Slah Drira, Ralf Frotscher, Manfred Staat. An element-based formulation for ES-FEM and FS-FEM models for implementation in standard solid mechanics finite element codes for 2D and 3D static analysis. International Journal for Numerical Methods in Engineering, 2023, 124 (2), pp.402-433. 10.1002/nme.7126. hal-03956203

HAL Id: hal-03956203

https://ifp.hal.science/hal-03956203

Submitted on 25 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

<figure><img src="image_2.png" /></figure>

# An element-based formulation for ES-FEM and FS-FEM models for implementation in standard solid mechanics finite element codes for 2D and 3D static analysis

Daniele Colombo<fnref n="1" /> ID | Slah Drira<fnref n="2" /> | Ralf Frotscher<fnref n="3" /> ID | Manfred Staat<fnref n="3" /> ID

**Correspondence**

Daniele Colombo, IFP Energies Nouvelles, 1-4 avenue de Bois-Préau, 92852 Rueil-Malmaison, France.
Email: daniele.colombo@ifpen.fr

**Abstract**

Edge-based and face-based smoothed finite element methods (ES-FEM and FS-FEM, respectively) are modified versions of the finite element method allowing to achieve more accurate results and to reduce sensitivity to mesh distortion, at least for linear elements. These properties make the two methods very attractive. However, their implementation in a standard finite element code is nontrivial because it requires heavy and extensive modifications to the code architecture. In this article, we present an element-based formulation of ES-FEM and FS-FEM methods allowing to implement the two methods in a standard finite element code with no modifications to its architecture. Moreover, the element-based formulation permits to easily manage any type of element, especially in 3D models where, to the best of the authors' knowledge, only tetrahedral elements are used in FS-FEM applications found in the literature. Shape functions for non-simplex 3D elements are proposed in order to apply FS-FEM to any standard finite element.

**KEYWORDS**

distorted element, ES-FEM, FS-FEM, non-simplex S-FEM elements, S-FEM, smoothed finite element method

1 | INTRODUCTION

The smoothed finite element method (S-FEM) is a modified version of the standard finite element method (FEM) characterized by interesting properties making it a very attractive method. For a given mesh, S-FEM allows to achieve more accurate results than FEM, is less sensitive to mesh distortion and can be easily applied to polygonal and polyhedral elements, since element shape function values have to be known at some particular points in an element and no analytical form is required.<fnref n="1" />⁻<fnref n="3" />

From a practical point of view, S-FEM is based on the smoothed strain field that is calculated by averaging standard FEM compatible strain field over so-called smoothing domains. These domains are obtained by subdividing an FEM mesh into nonoverlapping domains. Construction of smoothing domains can be accomplished in several ways. The features and properties of the resulting S-FEM method are strongly influenced by the adopted construction procedure. Four S-FEM

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2022 The Authors. *International Journal for Numerical Methods in Engineering* published by John Wiley & Sons Ltd.

<footnote n="1">¹Département Physique Numérique des Milieux Poreux, IFP Energies Nouvelles, Rueil-Malmaison, France</footnote>

<footnote n="2">²Laboratoire d'informatique et de mécanique appliquées à la construction - IMAC, Ecole Polytechnique Fédérale de Lausanne - EPFL, Lausanne, Switzerland</footnote>

<footnote n="3">³Faculty of Medical Engineering and Applied Mathematics, FH Aachen University of Applied Sciences, Jülich, Germany</footnote>
methods are currently available in the literature, each one using a different geometrical mesh entity to build the smoothing domains:¹,² CS-FEM (cell-based S-FEM, using the whole element),⁴ NS-FEM (nodal-based S-FEM, using mesh nodes),⁵ ES-FEM (edge-based S-FEM, using 2D elements edges)⁶⁻⁹ for 2D models only, and FS-FEM (face-based S-FEM, using 3D element faces) for 3D models only.¹⁰

Among these S-FEM methods, the ES-FEM (2D models) and the FS-FEM (3D models) are the most attractive because they achieve very accurate results and high mesh convergence rates compared to FEM.² However, their implementation in a standard finite element code is not straightforward. Discretization of solid mechanics equations by means of the S-FEM method results in a set of algebraic equations with the same form of the one obtained by the standard finite element procedure. The unknowns are still associated to mesh nodes and they are exactly the same as in the case of the standard finite element discretization. However, the stiffness matrix is no more calculated on an element basis but on a smoothing domain basis: a stiffness matrix is calculated for each smoothing domain and is assembled in the global stiffness matrix of the mesh by an assembly process similar to the one used in standard FEM.

The fact that in ES-FEM and FS-FEM methods a smoothing domain extends over two adjacent elements makes the implementation of S-FEM in a standard FEM code a challenging task, since extensive and substantial modifications of the code are required in order to account for the new entities that are used to calculate the global stiffness matrix contributions. In a standard FE code, element stiffness matrix calculations and their assembly in the global stiffness matrix of the model are done element by element. A new element formulation is added to the code by implementing all the calculations for one element only. In the ES-FEM and FS-FEM methods, stiffness matrix calculations are done in one smoothing domain, which extends over two adjacent elements. All the degrees of freedom of the two elements are used and information from the two adjacent elements are necessary. Therefore, it is impossible to implement an ES-FEM or FS-FEM element by simply adding a new element in a standard FE code because, due to the code architecture, this element has no access to information from the adjacent elements and only the degrees of freedom of the element are used for the calculation of the stiffness matrix.

In this article, an element-based formulation of ES-FEM and FS-FEM methods is presented. This formulation allows the implementation of ES-FEM and FS-FEM methods in a standard finite element code by using the standard procedure usually adopted to add a new standard FEM element. This is achieved by expressing the smoothing domain stiffness matrix as the sum of the contributions coming from two adjacent elements containing the smoothing domain and from an additional “auxiliary element” associated to each smoothing domain. Such a strategy permits to execute all the smoothing domain stiffness matrix calculations on an element basis, while avoiding any modification to the standard finite element code architecture.

The possibility to implement the element-based formulation in a standard FE code depends on the capabilities offered by the user interface of the code. For an open-source code, the user has access to any part of the code and she or he can freely implement the formulation. For a closed-source code, the user can successfully implement the formulation if the interface to the code allows to:

* Define additional elements starting from a given mesh in order to define the auxiliary elements used by the element-based formulation.

* Implement stiffness matrix and strain calculations for a user-defined element.

* Exchange information between the user defined elements in order to permit to one element to get data from the adjacent elements.

For example, the widely used Abaqus solver offers to the user powerful tools to extend code capabilities. The user could use Abaqus scripting to create auxiliary elements for a given input mesh, could define an Abaqus user element to implement the stiffness matrix calculation for the element-based formulation for the elements and auxiliary elements of the mesh and could use Abaqus scripting or Abaqus user subroutines to exchange data between the elements. Data used by the element-based formulation could be stored in state variables of the user-defined elements.

It is clear that if the possibility to implement the element-based S-FEM formulation in a closed-source standard FE code depends on the user interface capabilities of the code, the standard S-FEM formulation cannot be implemented because the user has no mean to modify the code architecture through the user interface of the code, independently of its capabilities. These modifications are necessary for the standard S-FEM formulation, as discussed above.
This article is organized as follows. In Section 2, the governing equations of the static solid mechanics problems are briefly reviewed. In Section 3, the basic equations produced by the standard finite element discretization are presented. In Section 4, these equations are compared to the equations produced by S-FEM discretization in order to highlight the ideas underlying the element-based formulation for ES-FEM and FS-FEM methods, which is presented in detail in Section 5 with the shape function formulation for non-simplex elements. Finally, the element-based formulation is applied to a 2D and a 3D case in Section 6 and the results are compared to FEM results and reference solutions available in literature.

## 2 | GOVERNING EQUATIONS

We consider a static solid mechanic problem defined in a domain $\Omega \in \mathbb{R}^m$ subjected to external body forces $\mathbf{b} = [b_1 \dots b_m]^T$. Prescribed displacements $\mathbf{u}_\Gamma = [u_{\Gamma_1} \dots u_{\Gamma_m}]^T$ and external surface forces $\mathbf{t}_\Gamma = [t_{\Gamma_1} \dots t_{\Gamma_m}]^T$ are applied on the displacement and traction boundary on $\Gamma_u \in \Gamma = \partial\Omega$ and $\Gamma_t \in \Gamma = \partial\Omega$, respectively.

The equilibrium equation can be written as:2,11,12

$$ \frac{\partial \sigma_{ij}}{\partial x_j} + b_i = 0, \quad i, j = 1, \dots, m \quad \text{for } x_j \in \Omega \setminus \partial \Omega, \qquad (1) $$

where $\sigma_{ij}$ are the stress tensor components and $\mathbf{x} = [x_1 \dots x_m]^T$ is the position vector defining point position. Equation (1) can be equivalently written in matrix form as follows:

$$ \mathbf{L}_d^T \boldsymbol{\sigma} + \mathbf{b} = \mathbf{0}, \qquad (2) $$

where $\mathbf{L}_d$ is the matrix of differential operators given by:

$$ \mathbf{L}_d = \begin{bmatrix} \frac{\partial}{\partial x_1} & 0 \\ 0 & \frac{\partial}{\partial x_2} \\ \frac{\partial}{\partial x_2} & \frac{\partial}{\partial x_1} \end{bmatrix} \quad \text{for 2D case (m = 2)}, \\ \mathbf{L}_d = \begin{bmatrix} \frac{\partial}{\partial x_1} & 0 & 0 \\ 0 & \frac{\partial}{\partial x_2} & 0 \\ 0 & 0 & \frac{\partial}{\partial x_3} \\ 0 & \frac{\partial}{\partial x_3} & \frac{\partial}{\partial x_2} \\ \frac{\partial}{\partial x_3} & 0 & \frac{\partial}{\partial x_1} \\ \frac{\partial}{\partial x_2} & \frac{\partial}{\partial x_1} & 0 \end{bmatrix} \quad \text{for 3D case (m = 3)}, \qquad (3) $$

and $\boldsymbol{\sigma}$ is the stress vector given by:

$$ \boldsymbol{\sigma} = [\sigma_{11} \ \sigma_{22} \ \sigma_{12}]^T \quad \text{for 2D case (m = 2)}, \\ \boldsymbol{\sigma} = [\sigma_{11} \ \sigma_{22} \ \sigma_{33} \ \sigma_{23} \ \sigma_{13} \ \sigma_{12}]^T \quad \text{for 3D case (m = 3)}. \qquad (4) $$

The strain field is defined by means of the compatibility equation under the assumption of small displacements:2,11,12

$$ \epsilon_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right), \quad i, j = 1, \dots, m \quad \text{for } x_j \in \Omega \setminus \partial \Omega, \qquad (5) $$

where $u_i$ are the components of the displacement vector $\mathbf{u} = [u_1 \dots u_m]^T$. Equation (5) can be equivalently written in matrix form as:

$$ \boldsymbol{\epsilon} = \mathbf{L}_d \mathbf{u}, \qquad (6) $$
where $\epsilon$ is the strain vector given by:

$$ 
\begin{aligned}
\epsilon &= [\epsilon_{11} \ \epsilon_{22} \ 2\epsilon_{12}]^T = [\epsilon_{11} \ \epsilon_{22} \ \gamma_{12}]^T \quad \text{for 2D case (m = 2)}, \\
\epsilon &= [\epsilon_{11} \ \epsilon_{22} \ \epsilon_{33} \ 2\epsilon_{23} \ 2\epsilon_{13} \ 2\epsilon_{12}]^T = [\epsilon_{11} \ \epsilon_{22} \ \epsilon_{33} \ \gamma_{23} \ \gamma_{13} \ \gamma_{12}]^T \quad \text{for 3D case (m = 3)}.
\end{aligned}
\quad (7) $$

The constitutive equation relating stress and strain is given by the generalized Hooke's law:²,¹¹,¹²

$$ \sigma_{ij} = c_{ijkl}\epsilon_{kl}, \quad i, j, k, l \in [1, m] \quad (8) $$

or in equivalent matrix form:

$$ \sigma = C\epsilon. \quad (9) $$

The fourth-order tensor $C_{ijkl}$ refers to material stiffness moduli.

Finally, the boundary conditions read:

$$ 
\begin{aligned}
\mathbf{u} &= \mathbf{u}_{\Gamma_u} \quad \text{on } \Gamma_u, \\
\mathbf{L}_n^T \boldsymbol{\sigma} &= \mathbf{t}_{\Gamma_t} \quad \text{on } \Gamma_t,
\end{aligned}
\quad (10a) $$

where $\mathbf{L}_n(\mathbf{x})$ is the matrix of the outward normals $\mathbf{n}(\mathbf{x})$ to the domain boundary $\partial\Omega$ given by:

$$ 
\begin{aligned}
\mathbf{L}_n(\mathbf{x}) &= \begin{bmatrix} n_x(\mathbf{x}) & 0 \\ 0 & n_y(\mathbf{x}) \\ n_y(\mathbf{x}) & n_x(\mathbf{x}) \end{bmatrix} && \text{for 2D case, } \mathbf{n}(\mathbf{x}) = [n_x(\mathbf{x}) \ n_y(\mathbf{x})]^T, \\
\mathbf{L}_n(\mathbf{x}) &= \begin{bmatrix} n_x(\mathbf{x}) & 0 & 0 \\ 0 & n_y(\mathbf{x}) & 0 \\ 0 & 0 & n_z(\mathbf{x}) \\ 0 & n_z(\mathbf{x}) & n_y(\mathbf{x}) \\ n_z(\mathbf{x}) & 0 & n_x(\mathbf{x}) \\ n_y(\mathbf{x}) & n_x(\mathbf{x}) & 0 \end{bmatrix} && \text{for 3D case, } \mathbf{n}(\mathbf{x}) = [n_x(\mathbf{x}) \ n_y(\mathbf{x}) \ n_z(\mathbf{x})]^T.
\end{aligned}
\quad (10b) $$

### 3 | FEM DISCRETIZATION

A finite element mesh is built on the domain $\Omega$. The mesh is composed by $n_e$ elements and $n_n$ nodes. Displacement degrees of freedom are assigned to each node $i$:

$$ \mathbf{d}_i = [d_1 \ \dots \ d_m]^T. \quad (11) $$

The global vector of the degrees of freedom of the whole mesh is given by:

$$ \mathbf{d} = [\mathbf{d}_1^T \ \dots \ \mathbf{d}_{n_n}^T]^T. \quad (12) $$

The displacement field is approximated by means of element shape functions:¹¹⁻¹³

$$ \mathbf{u}(\mathbf{x}) = \mathbf{N}(\mathbf{x})\mathbf{d}, \quad (13) $$

where $\mathbf{N}(\mathbf{x}) = [\mathbf{N}_1 \ \dots \ \mathbf{N}_{n_n}]$ is the global shape function matrix, $\mathbf{N}_i(\mathbf{x}) = N_i(\mathbf{x}) \mathbf{I}_m$ is the shape function matrix for node $i$, $N_i(\mathbf{x})$ is the shape function for the same node, and $\mathbf{I}_m$ is the identity matrix of size $m$. As a consequence of its definition, matrix $\mathbf{N}(\mathbf{x})$ is coarse.
The strain field is therefore approximated as follows (Equations 6 and 13):¹¹⁻¹³

$$ \varepsilon = L_d N(\mathbf{x}) d = B(\mathbf{x}) d, \qquad (14) $$

where $B(\mathbf{x}) = L_d N(\mathbf{x})$ is the strain-displacement matrix. For the following discussion, it is worth rewriting $B(\mathbf{x})$ in a form where contributions of each node appear explicitly:

$$ B(\mathbf{x}) = [B_1(\mathbf{x}) \dots B_{n_n}(\mathbf{x})], \qquad (15) $$

where $B_i(\mathbf{x}) = L_d N_i(\mathbf{x})$ is the strain-displacement matrix for node i.

FE method solves equilibrium equation (1) in weak form adopting the above discretized displacement and strain fields (Equations 13 and 14). The resulting system of algebraic equations can be written in matrix form as follow:¹¹⁻¹³

$$ Kd = f, \qquad (16) $$

where $K$ is the global stiffness matrix given by:

$$ K = \int_{\Omega} B^{T}(\mathbf{x}) CB(\mathbf{x}) d\Omega, \qquad (17) $$

and $f$ is the vector of nodal forces given by:

$$ f = \int_{\Omega} N^{T}(\mathbf{x}) b d\Omega + \int_{\Gamma_{i}} N^{T}(\mathbf{x}) t_{\Gamma} d\Gamma. \qquad (18) $$

For the following discussion, it is worth noting that contributions of each finite element to the global stiffness matrix can be easily separated. Indeed the global stiffness matrix is very coarse, which, as already mentioned, is due to the fact that matrix $N(\mathbf{x})$ is coarse. Equation (17) can be easily rewritten on an element-basis:

$$ K = \sum_{i=1}^{n_e} \int_{\Omega_i} B^{T}(\mathbf{x}) CB(\mathbf{x}) d\Omega. \qquad (19) $$

Considering the fact that integration is now carried out on a single element $\Omega_i$ and that only the shape functions of the element nodes contribute to the integrand, the global stiffness matrix can be equivalently rewritten as the sum of the stiffness matrices of each element $\Omega_i$ (“summation” here means “assembly” since the involved matrices have different size, see below for details):

$$ K = \sum_{i=1}^{n_e} K_{\Omega_i}, \qquad (20) $$

where the element stiffness matrix is defined as:

$$ K_{\Omega_i} = \int_{\Omega_i} B_{\Omega_i}^{T}(\mathbf{x}) CB_{\Omega_i}(\mathbf{x}) d\Omega. \qquad (21) $$

Strain-displacement matrix $B_{\Omega_i}(\mathbf{x})$ is the restriction to element $\Omega_i$ of the global strain-displacement matrix $B(\mathbf{x})$ (Equation 15) and is given by:

$$ B_{\Omega_i}(\mathbf{x}) = [B_1(\mathbf{x}) \dots B_{n_{\Omega_i}}(\mathbf{x})], \qquad (22) $$

where $n_{\Omega_i}$ is the number of nodes of element $\Omega_i$.

The size of element stiffness matrix $K_{\Omega_i}$ is $(mn_{\Omega_i}) \times (mn_{\Omega_i})$, where $n_{\Omega_i}$ is the number of nodes of element $\Omega_i$ and $m$ is the dimension of domain $\Omega$. On the contrary, the size of the global stiffness matrix $K$ is $(mn_n) \times (mn_n)$, where $n_n$ is the
number of nodes of the model. Strictly speaking, the summation symbol in Equation (20) is used improperly, because the size of the matrices to sum is not the same. Equation (20) is just a convenient way to concisely express the assembly procedure commonly used in finite element codes: the element stiffness matrix is first calculated for an element and then added to the global stiffness matrix by selecting the rows and columns corresponding to the degrees of freedom of the finite element.¹²,¹³

The element-based assembly procedure will be exploited later to define the S-FEM implementation strategy.

## 4 | S-FEM DISCRETIZATION

S-FEM models work on the so-called smoothing domains obtained by subdividing the finite element mesh of domain $\Omega$ into several nonoverlapping and no gap subdomains. A domain can be contained entirely inside one element or can extend over two adjacent elements, depending on the type of S-FEM method. For each smoothing domain $\Omega_k$, a smoothed strain is calculated at any point $\mathbf{x}_k \in \Omega_k$ by means of the strain/gradient smoothing technique:²

$$ \varepsilon_{\Omega_k}(\mathbf{x}_k) = \int_{\Omega_k} \mathbf{L}_d(\mathbf{x}) \mathbf{u}(\mathbf{x}) W(\mathbf{x}_k - \mathbf{x}) d\Omega, \quad (23) $$

where $W(\mathbf{x}_k - \mathbf{x})$ is a smoothing function defined as:

$$ W(\mathbf{x}_k - \mathbf{x}) = \begin{cases} \frac{1}{S_{\Omega_k}}, & \text{for } \mathbf{x} \in \Omega_k, \\ 0, & \text{for } \mathbf{x} \notin \Omega_k, \end{cases} \quad (24) $$

where $S_{\Omega_k}$ is the geometrical size of domain $\Omega_k$ (its area in 2D and its volume in 3D).

Because the displacement field $\mathbf{u}(\mathbf{x})$ is assumed continuous on $\partial\Omega_k$ and the smoothing function $W(\mathbf{x}_k - \mathbf{x})$ is differentiable over $\Omega_k$, by applying Green's divergence theorem, the smoothed strain can be rewritten equivalently as:²

$$ \varepsilon_{\Omega_k}(\mathbf{x}_k) = \varepsilon_{\Omega_k} = \frac{1}{S_{\Omega_k}} \int_{\partial\Omega_k} \mathbf{L}_n(\mathbf{x}) \mathbf{u}(\mathbf{x}) d\Gamma. \quad (25) $$

In S-FEM, it is assumed that the strain $\varepsilon_{\Omega_k}$ in the smoothing domain $\Omega_k$ is constant and equal to $\varepsilon_{\Omega_k}(\mathbf{x}_k)$ given by Equation (25).²,⁶,¹⁰,¹⁴

The integral in Equation (24) is restricted to the smoothing domain boundary $\partial\Omega_k$.

The discretization of the displacement field is identical to the one used in the FE method (Equation 13). The smoothed strain (Equation 24) can therefore be rewritten in matrix form as:²,⁶,¹⁰,¹⁴

$$ \varepsilon_{\Omega_k} = \mathbf{B}_{\Omega_k} \mathbf{d}, \quad (26) $$

where $\mathbf{B}_{\Omega_k}$ is the strain-displacement matrix associated to smoothing domain $\Omega_k$, given by:

$$ \mathbf{B}_{\Omega_k} = [\mathbf{B}_1^{\Omega_k} \dots \mathbf{B}_r^{\Omega_k}], \quad (27) $$

where $r$ is the total number of nodes of the elements containing the smoothing domain. The contribution of each node is given by $\mathbf{B}_i^{\Omega_k}$, which is the strain-displacement matrix of node $i$ associated to smoothing domain $\Omega_k$, given by:

$$ \mathbf{B}_i^{\Omega_k} = \frac{1}{S_{\Omega_k}} \int_{\partial\Omega_k} \mathbf{L}_n(\mathbf{x}) \mathbf{N}_i(\mathbf{x}) d\Gamma. \quad (28) $$

Matrix $\mathbf{B}_i^{\Omega_k}$ has the same structure as matrix $\mathbf{N}_i(\mathbf{x})$:²,¹⁰

$$ \mathbf{B}_i^{\Omega_k} = \begin{bmatrix} b_{i1}^{\Omega_k} & 0 \\ 0 & b_{i2}^{\Omega_k} \\ b_{i2}^{\Omega_k} & b_{i1}^{\Omega_k} \end{bmatrix} \quad \text{for 2D case,} \quad \mathbf{b}_i^{\Omega_k} = \begin{bmatrix} b_{i1}^{\Omega_k} & b_{i2}^{\Omega_k} \end{bmatrix}^T, $$
$$ \mathbf{B}_i^{\Omega_k} = \begin{bmatrix} b_{i1}^{\Omega_k} & 0 & 0 \\ 0 & b_{i2}^{\Omega_k} & 0 \\ 0 & 0 & b_{i3}^{\Omega_k} \\ 0 & b_{i3}^{\Omega_k} & b_{i2}^{\Omega_k} \\ b_{i3}^{\Omega_k} & 0 & b_{i1}^{\Omega_k} \\ b_{i2}^{\Omega_k} & b_{i1}^{\Omega_k} & 0 \end{bmatrix} \quad \text{for 3D case,} \quad \mathbf{b}_i^{\Omega_k} = [b_{i1}^{\Omega_k} \ b_{i2}^{\Omega_k} \ b_{i3}^{\Omega_k}]^T, \qquad (29) $$

where

$$ b_{ij}^{\Omega_k} = \frac{1}{S_{\Omega_k}} \int_{\partial\Omega_k} n_j(\mathbf{x}) N_i(\mathbf{x}) d\Gamma. \qquad (30) $$

If linear elements are used, the shape functions are linear and the integration in Equation (30) is trivial:²,¹⁰

$$ b_{ij}^{\Omega_k} = \frac{1}{S_{\Omega_k}} \sum_{p=1}^{n_{\Gamma}} n_j(\mathbf{x}_{\text{Gp}}) N_i(\mathbf{x}_{\text{Gp}}) l_p, \qquad (31) $$

where $n_Γ$ is the number of segments (2D case) or planar faces (3D case) constituting the smoothing domain boundary $\partial\Omega_k$, $\mathbf{x}_{\text{Gp}}$ is the position of the $p$th segment mid-point (2D case) or the $p$th planar face centroid (3D case), and $l_p$ is its length (2D case) or its area (3D case).

It can be demonstrated that a smoothed Galerkin weak form of the equilibrium equation (1) based on smoothed strain can be written as done in FE case, leading to the following algebraic equation system:²

$$ \mathbf{K}_{\text{SFEM}} \mathbf{d} = \mathbf{f}. \qquad (32) $$

Equation (32) is formally identical to Equation (16) obtained for FEM discretization. Vector **f** is given by Equation (18) as in the FEM case. The global stiffness matrix **K**<sub>SFEM</sub> is formally obtained by exactly the same assembly procedure as used in FEM but based on the contributions of each smoothing domain (Equation 20):

$$ \mathbf{K}_{\text{SFEM}} = \sum_{k=1}^{n_{\Omega}} \mathbf{K}_{\Omega_k \text{SFEM}}, \qquad (33) $$

where $n_Ω$ is the number of smoothing domains. The smoothing domain stiffness matrix $\mathbf{K}_{\Omega_k \text{SFEM}}$ is given by the following equation, which is formally identical to Equation (17) used in the FEM case:

$$ \mathbf{K}_{\Omega_k \text{SFEM}} = \int_{\Omega_k} \mathbf{B}_{\Omega_k}^T \mathbf{C} \mathbf{B}_{\Omega_k} d\Omega = \mathbf{B}_{\Omega_k}^T \mathbf{C} \mathbf{B}_{\Omega_k} S_{\Omega_k}. \qquad (34) $$

Equation (34) shows that the stiffness matrix calculation does not need integration and it can be computed directly from the values of $\mathbf{B}_{\Omega_k} = [\mathbf{B}_1^{\Omega_k} \dots \mathbf{B}_r^{\Omega_k}]$ and $\mathbf{C}$ matrices. Equation (31) shows that, if linear shape functions are used on the smoothing domain boundary $\partial\Omega_k$, integration of terms in matrix $\mathbf{B}_i^{\Omega_k}$ are done directly in the physical space and no isoparametric transformation and no evaluation of the Jacobian of the associated transformation matrix are needed, contrary to FEM (see Equation 21). This makes S-FEM method more insensitive to mesh distortion,² as it will be shown in Section 6. Formulation of linear shape functions on the smoothing domain boundary $\partial\Omega_k$ for non-simplex elements will be done in Section 5.

Extension of the S-FEM formulation to large deformations and geometrically nonlinear problems can be done with a similar procedure to the one used in standard FEM, as detailed in References 2 and 10. A straightforward application of S-FEM formulation to this class of problems can be done by adopting an updated Lagrangian approach.¹¹ Discussion of this topic is out of the scope of this article and interested readers will find all the details in References 2, 10, and 11.
# 5 | AN ELEMENT-BASED FORMULATION FOR ES-FEM AND FS-FEM METHODS

Equations (20) and (33) are formally identical, except that the assembly procedure is carried out on smoothing domains for S-FEM, whereas it is carried out on finite elements for FEM. The architecture of finite element codes is usually based on elements and the calculations for the element stiffness matrices $K_{\Omega_i}$ (Equation 21) are carried out element by element without the need of exchanging data between elements. Thus the parallelization of these calculations is straightforward.

The same approach can be used for S-FEM. However, all S-FEM methods except CS-FEM create smoothing domains that extend over two adjacent elements (Figures 1 and 2 for ES-FEM and FS-FEM methods, respectively) and data from the two elements are needed to calculate the smoothing domain stiffness matrix $K_{\Omega_k}$ in Equation (33). Consequently, contrary to the FEM, the degrees of freedom of all the nodes of the two adjacent elements are related in $K_{\Omega_k}$, Therefore the implementation of NS-FEM, ES-FEM and FS-FEM methods in an existing finite element code by means of Equation (33) requires modifications to the code architecture since FEM codes usually adopt an element-based formulation. These modifications are very complex and time-consuming.

In the following paragraphs, it will be shown that S-FEM equations can be equivalently rewritten using an element-based formulation for ES-FEM and FS-FEM methods. This equivalent formulation allows easy implementation of the S-FEM method into an existing FEM code because it does not require modifications of code architecture.

## 5.1 | Element-based formulation for smoothed stiffness matrix $K_{\Omega_k}$ for ES-FEM and FS-FEM methods

Equation (33) can be rewritten in a different way in order to identify the contribution of each element containing the smoothing domain.

We consider a smoothing domain $\Omega_k$ extending over two adjacent elements A and B. Some nodes are shared between these two elements. In 2D, the ES-FEM method builds the smoothing domain at the common edge between the two elements. Consequently, the two nodes defining the edge are shared between element A and B. In 3D, FS-FEM method builds the smoothing domain on the common face between the two elements. Therefore, for commonly used finite elements,

<figure><img src="image_3.png" /><figcaption>**FIGURE 1** Example of construction of a smoothing domain for ES-FEM method. The smoothing domain is built on the edge defined by nodes 2 and 4 shared between the quadrangular and triangular elements by connecting these two nodes with the element centroids A and B.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**FIGURE 2** Example of construction of a smoothing domain for FS-FEM method. The smoothing domain is built on the face defined by nodes 5, 6, 7, and 8 shared between the two hexahedral elements by connecting these four nodes with the element centroids A and B.</figcaption></figure>
three (triangular face) or four (quadrangular face) nodes are shared between element A and B. For polyhedral elements, a greater number of nodes is shared.

By rearranging the node order, contributions to the strain-displacement matrix $B_{\Omega_k}$ of the nodes belonging to one element only and of the shared nodes can be separated in Equation (27):

$$B_{\Omega_k} = \left[ \begin{array}{c|c|c} \text{element A only} & \text{shared A and B} & \text{element B only} \\ \mathbf{B}_{1}^{\Omega_k} \dots \mathbf{B}_{n_A}^{\Omega_k} | \mathbf{B}_{n_A+1}^{\Omega_k} \dots \mathbf{B}_{n_A+n_S}^{\Omega_k} | \mathbf{B}_{n_A+n_S+1}^{\Omega_k} \dots \mathbf{B}_{n_A+n_S+n_B}^{\Omega_k} \end{array} \right], \quad (35)$$

where $n_A$ and $n_B$ are the number of nodes belonging only to element A and B, respectively, and $n_S$ is the number of shared nodes between the two elements.

We denote $B_{\Omega_k}|_A$, $B_{\Omega_k}|_B$, and $B_{\Omega_k}|_S$ the smoothing domain strain-displacement matrix $B_{\Omega_k}$ restricted to the nodes belonging only to elements A and B and shared between element A and B, respectively:

$$ \begin{aligned} B_{\Omega_k}|_A &= [\mathbf{B}_{1}^{\Omega_k} \dots \mathbf{B}_{n_A}^{\Omega_k}], \\ B_{\Omega_k}|_S &= [\mathbf{B}_{n_A+1}^{\Omega_k} \dots \mathbf{B}_{n_A+n_S}^{\Omega_k}], \\ B_{\Omega_k}|_B &= [\mathbf{B}_{n_A+n_S+1}^{\Omega_k} \dots \mathbf{B}_{n_A+n_S+n_B}^{\Omega_k}]. \end{aligned} \quad (36) $$

By substituting Equation (36) in Equation (35), matrix $B_{\Omega_k}$ can be equivalently written as:

$$B_{\Omega_k} = [B_{\Omega_k}|_A B_{\Omega_k}|_S B_{\Omega_k}|_B]. \quad (37)$$

By substituting Equation (37) into Equation (34), the smoothing domain stiffness matrix $K_{\Omega_k}$ can be written as:

$$K_{\Omega_k} = [B_{\Omega_k}|_A B_{\Omega_k}|_S B_{\Omega_k}|_B]^T C [B_{\Omega_k}|_A B_{\Omega_k}|_S B_{\Omega_k}|_B] S_{\Omega_k} = \begin{bmatrix} K_{AA}^{\Omega_k} & K_{AS}^{\Omega_k} & K_{AB}^{\Omega_k} \\ K_{SA}^{\Omega_k} & K_{SS}^{\Omega_k} & K_{SB}^{\Omega_k} \\ K_{BA}^{\Omega_k} & K_{BS}^{\Omega_k} & K_{BB}^{\Omega_k} \end{bmatrix} S_{\Omega_k}, \quad (38)$$

where $K_{IJ}^{\Omega_k}$ is defined as:

$$K_{IJ}^{\Omega_k} = (B_{\Omega_k}|_I)^T C B_{\Omega_k}|_J, \quad \text{for } I, J \in \{A, B, S\}. \quad (39)$$

Contributions of each element to $K_{\Omega_k}$ in Equation (38) can now be easily separated:

$$K_{\Omega_k} = \left[ \begin{array}{c|c} \text{element A} & \text{element B} \\ \begin{bmatrix} \mathbf{K}_{AA}^{\Omega_k} & \mathbf{K}_{AS}^{\Omega_k} \\ \mathbf{K}_{SA}^{\Omega_k} & \frac{\mathbf{K}_{SS}^{\Omega_k}}{2} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} & \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \frac{\mathbf{K}_{SS}^{\Omega_k}}{2} & \mathbf{K}_{SB}^{\Omega_k} \\ \mathbf{0} & \mathbf{K}_{BS}^{\Omega_k} & \mathbf{K}_{BB}^{\Omega_k} \end{bmatrix} + \begin{bmatrix} \mathbf{0} & \mathbf{0} & \mathbf{K}_{AB}^{\Omega_k} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{K}_{BA}^{\Omega_k} & \mathbf{0} & \mathbf{0} \end{bmatrix} \end{array} \right]. \quad (40)$$

In Equation (40), $K_{\Omega_k}$ is expressed as the sum of three matrices obtained by grouping the terms that can be calculated by accounting for the nodes of one element only: element A can calculate the terms $K_{AA}^{\Omega_k}$, $K_{AS}^{\Omega_k}$, $K_{SA}^{\Omega_k}$, and $K_{SS}^{\Omega_k}$ involving the nodes of element A only and the shared nodes, and similarly element B can calculate the terms $K_{BB}^{\Omega_k}$, $K_{BS}^{\Omega_k}$, $K_{SB}^{\Omega_k}$, and $K_{SS}^{\Omega_k}$. Equation (40) clearly shows that the two terms $K_{AB}^{\Omega_k}$ and $K_{BA}^{\Omega_k}$ cannot be calculated by considering the two elements separately: these terms relate nodes belonging to element A and B only. In order to calculate these terms, a new element must be added to the mesh. We will call this element the auxiliary element of smoothing domain $\Omega_k$. A detailed discussion of its properties will be done later in Section 5.3.

Submatrix $K_{SS}^{\Omega_k}$ in Equation (38) needs a special treatment. Its elements are calculated from the contributions to the shared nodes of both elements A and B. This involves data exchange between the two elements because each element can only do a partial calculation, as it will be detailed in Section 5.2.1. Both elements A and B can evaluate $K_{SS}^{\Omega_k}$. Therefore the
calculated value must be divided by two, since the contributions of elements A and B are added to obtain the final stiffness matrix $K_{\Omega_k}$, as shown in Equation (40). Alternatively, a more effective solution, which can be implemented without additional computational costs, would be to calculate $K_{\Omega_k}^{\Omega_k}$ in only one element, avoiding to make the same calculations twice. Any of the two elements A and B can do the calculation. The choice is arbitrary and can be done during the creation of the smoothing domains at the beginning of the simulation. For example, if element A has been chosen for the calculation of $K_{\Omega_k}^{\Omega_k}$ submatrix, Equation (40) would be rewritten as:

$$ K_{\Omega_k} = \left[ \begin{array}{c|c|c} \text{element A} & \text{element B} & \text{auxiliary element} \\ \begin{bmatrix} K_{AA}^{\Omega_k} & K_{AS}^{\Omega_k} & 0 \\ K_{SA}^{\Omega_k} & K_{SS}^{\Omega_k} & 0 \\ 0 & 0 & 0 \end{bmatrix} & \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & K_{SB}^{\Omega_k} \\ 0 & K_{BS}^{\Omega_k} & K_{BB}^{\Omega_k} \end{bmatrix} & \begin{bmatrix} 0 & 0 & K_{AB}^{\Omega_k} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{array} \right]. \quad (41) $$

## 5.2 | Evaluation of $K_{IJ}^{\Omega_k}$ submatrices

The stiffness matrices of elements A and B and auxiliary element (Equation 41) are built using submatrices $K_{IJ}^{\Omega_k}$ given by Equation (39). These submatrices are calculated from strain-displacement matrices $B_{\Omega_k}|_A$, $B_{\Omega_k}|_B$, and $B_{\Omega_k}|_S$, which are a restriction of strain-displacement matrix $B_{\Omega_k}$ to the nodes belonging to elements A and B exclusively and to the shared nodes, respectively (Equation 36). Matrices $B_{\Omega_k}|_A$, $B_{\Omega_k}|_B$, and $B_{\Omega_k}|_S$ are therefore built using the strain-displacement matrices $B_i^{\Omega_k}$ of the nodes of the elements containing the smoothing domain $\Omega_k$ (Equation 28).

The elements of $B_i^{\Omega_k}$ are calculated by means of the shape function of node i on the boundary of the smoothing domain (Equations 28,30, and 31). As a consequence:

* If node i belongs exclusively to element E (element A or element B), $B_i^{\Omega_k}$ can be calculated by considering that element only. This is a consequence of the fact that shape function $N_i(x)$ is zero inside the element adjacent to element E and therefore its contributions to $B_i^{\Omega_k}$ are zero. Matrices $B_{\Omega_k}|_A$ and $B_{\Omega_k}|_B$ can be calculated considering only element A and B, respectively, that is no information from the adjacent element are necessary.

* If node i is a shared node, both elements A and B contribute to $B_i^{\Omega_k}$ because the shape function $N_i(x)$ is nonzero everywhere on the smoothing domain boundary. Therefore, the evaluation of $B_{\Omega_k}|_S$ needs information coming from the two elements A and B.

The element-based approach used in Equations (40) and (41) can therefore be only straightforwardly applied to the evaluation of $B_{\Omega_k}|_A$ and $B_{\Omega_k}|_B$ and further work is needed for $B_{\Omega_k}|_S$.

### 5.2.1 | Evaluation of $B_{\Omega_k}|_S$

Let us consider a node i that is shared between elements A and B. This node contributes to $B_{\Omega_k}|_S$ through the matrix $B_i^{\Omega_k}$ (Equation 36). Contributions of element A and B to $B_i^{\Omega_k}$ can be easily separated (Equation 28):

$$ B_i^{\Omega_k} = \frac{1}{S_{\Omega_k|A} + S_{\Omega_k|B}} \left( \underbrace{\int_{\partial\Omega_k|A} \mathbf{L}_n(\mathbf{x}) \mathbf{N}_i(\mathbf{x}) d\Gamma}_{\text{element A}} + \underbrace{\int_{\partial\Omega_k|B} \mathbf{L}_n(\mathbf{x}) \mathbf{N}_i(\mathbf{x}) d\Gamma}_{\text{element B}} \right) = \frac{1}{S_{\Omega_k|A} + S_{\Omega_k|B}} (B_i^{\Omega_k}|_A + B_i^{\Omega_k}|_B), \quad (42) $$

where $\Omega_k|A$ and $\Omega_k|B$ denote the restriction of the smoothing domain $\Omega_k$ to elements A and B, respectively. $\partial\Omega_k|A$ and $\partial\Omega_k|B$ denote the part of the boundary of smoothing domain $\Omega_k$ that is contained inside elements A and B, respectively, and $B_i^{\Omega_k}|_A$ and $B_i^{\Omega_k}|_B$ are defined as:

$$ B_i^{\Omega_k}|_E = \int_{\partial\Omega_k|E} \mathbf{L}_n(\mathbf{x}) \mathbf{N}_i(\mathbf{x}) d\Gamma \quad \text{for } E \in \{A, B\}. \quad (43) $$
Equation (42) clearly shows that element-based calculation of $B_i^{\Omega_k}$ must be carried out in three steps:

* First, each element E calculates $S_{\Omega_k|E}$ and $B_i^{\Omega_k}|E$. No information is needed to be exchanged between the two adjacent elements. This calculation can be straightforwardly parallelized by exploiting existing FEM code capabilities, since it is carried out on an element by element basis.

* Information is exchanged between the two adjacent elements: each element E sends $S_{\Omega_k|E}$ and $B_i^{\Omega_k}|E$ to the adjacent element. This step cannot be parallelized by exploiting the element-based parallelization of the FEM code.

* Finally, each element E calculates the value of $B_i^{\Omega_k}$ by means of Equation (42). This step can be parallelized by exploiting existing FEM code capabilities, since it is carried out on an element by element basis.

## 5.3 | Auxiliary element of smoothing domain $\Omega_k$

Equation (40) shows that an auxiliary element must be defined for each smoothing domain in order to correctly calculate its stiffness matrix.

The auxiliary element is built using, in the order, the nodes belonging to element A only and the nodes belonging to element B only, without using any of the shared nodes. It is only used to complete the stiffness matrix calculation by providing the stiffness matrices $K_{AB}^{\Omega_k}$ and $K_{BA}^{\Omega_k}$ (Equation 40). These matrices can be promptly built from the values of $B_i^{\Omega_k}$ coming from previous calculations done by elements A and B when evaluating $K_{AA}^{\Omega_k}$, $K_{AS}^{\Omega_k}$, $K_{SA}^{\Omega_k}$, $K_{BB}^{\Omega_k}$, $K_{BS}^{\Omega_k}$, $K_{SB}^{\Omega_k}$, and $K_{SS}^{\Omega_k}$.

* Matrix C and values of the vectors $b_i^{\Omega_k}$ are copied from elements A and B (Equations 29 and 30).

* The values $b_i^{\Omega_k}$ are used to build the matrices $B_i^{\Omega_k}$ (Equation 29) and then the matrices $B_{\Omega_k}|_A$ and $B_{\Omega_k}|_B$ (Equation 36), which are used together with matrix C to build the matrices $K_{AB}^{\Omega_k}$ and $K_{BA}^{\Omega_k}$ (equations 36 and 39) used to finally build the auxiliary element stiffness matrix. This step can be straightforwardly parallelized by exploiting existing FEM code capabilities, since it is carried out on an element by element basis.

The geometry of the auxiliary element is never used in the calculations, which makes its construction straightforward for any type of adjacent elements A and B. It is important to point out the fact that auxiliary elements are not used for any other calculation. They are not used in post-processing operations and they can therefore be omitted when displaying simulation results.

## 5.4 | Calculation of strain tensor for smoothing domain $\Omega_k$

Strain tensor is assumed to be constant over a smoothing domain $\Omega_k$ (Section 4). Its value is given by Equation (26), which can be rewritten as follows:

$$ \epsilon_{\Omega_k} = [B_{\Omega_k}|_A B_{\Omega_k}|_S B_{\Omega_k}|_B] [d_A^T d_S^T d_B^T]^T, \quad (44) $$

where the vectors $d_A$, $d_S$, and $d_B$ are the vectors containing the degrees of freedom of the nodes belonging to element A and element B exclusively and to shared nodes, respectively. They are defined as:

$$ \begin{aligned} d_A &= [d_1^T \dots d_{n_A}^T], \\ d_S &= [d_{n_A+1}^T \dots d_{n_A+n_S}^T], \\ d_B &= [d_{n_A+n_S+1}^T \dots d_{n_A+n_S+n_B}^T]. \end{aligned} \quad (45) $$

The contributions of each element to the strain tensor can be easily separated in Equation (44):

$$ \epsilon_{\Omega_k} = \underbrace{[B_{\Omega_k}|_A \frac{1}{2} B_{\Omega_k}|_S] [d_A^T d_S^T]^T}_{\text{element A}} + \underbrace{[\frac{1}{2} B_{\Omega_k}|_S B_{\Omega_k}|_B] [d_S^T d_B^T]^T}_{\text{element B}} \quad (46) $$

or in a more effective way, as done for stiffness matrix calculations (see Equation 41):
$$ \epsilon_{\Omega_k} = \overbrace{[\mathbf{B}_{\Omega_k}|_A \mathbf{B}_{\Omega_k}|_S]}^{\text{element A}} [\mathbf{d}_A^T \mathbf{d}_S^T]^T + \overbrace{[\mathbf{0} \ \mathbf{B}_{\Omega_k}|_B]}^{\text{element B}} [\mathbf{d}_S^T \mathbf{d}_B^T]^T. \quad (47) $$

As already pointed out in Section 5.3, auxiliary elements are not used for strain calculations.

Strain-displacement matrices $\mathbf{B}_{\Omega_k}|_A$, $\mathbf{B}_{\Omega_k}|_S$, and $\mathbf{B}_{\Omega_k}|_B$ have already been evaluated during stiffness matrix calculations (Section 5.2). In particular, matrix $\mathbf{B}_{\Omega_k}|_S$ for shared nodes has already been built and its value has already been stored for both elements A and B and therefore no further data exchange is needed between these elements.

Equation (47) clearly shows that the evaluation of the strain tensor has to be done in two steps:

* Elements A and B calculate their contribution to the smoothed strain tensor. This step can be straightforwardly parallelized by exploiting existing FEM code capabilities, since it is carried out on an element by element basis.

* The contributions of each element is exchanged with the adjacent element. Each element finally calculates the smoothed strain tensor by adding its contribution to the one coming from the adjacent element.

## 5.5 | Computation of internal force vector

In finite element discretization, the internal force vector $\mathbf{f}_{\text{int}}$ for element $\Omega_i$ is calculated as:¹¹,¹²

$$ \mathbf{f}_{\text{int}} = \int_{\Omega_i} \mathbf{B}_{\Omega_i}^T \sigma_{\Omega_i} d\Omega, \quad (48) $$

where $\mathbf{B}_{\Omega_i}$ is given by Equation (22).

In S-FEM discretization, the same equation holds for a smoothing domain $\Omega_k$:

$$ \mathbf{f}_{\text{int}} = \int_{\Omega_k} \mathbf{B}_{\Omega_k}^T \sigma_{\Omega_k} d\Omega = \mathbf{B}_{\Omega_k}^T \sigma_{\Omega_k} S_{\Omega_k}, \quad (49) $$

where $\mathbf{B}_{\Omega_k}$ is given by Equation (37). Stress tensor $\sigma_{\Omega_k}$ is constant over the smoothing domain $\Omega_k$.

By substituting Equation (37), the internal force vector can be equivalently rewritten as:

$$ \mathbf{f}_{\text{int}} = [\mathbf{B}_{\Omega_k}|_A \mathbf{B}_{\Omega_k}|_S \mathbf{B}_{\Omega_k}|_B]^T \sigma_{\Omega_k} S_{\Omega_k}. \quad (50) $$

Contributions of each element can now be easily separated:

$$ \mathbf{f}_{\text{int}} = \overbrace{[\mathbf{B}_{\Omega_k}|_A \frac{1}{2}\mathbf{B}_{\Omega_k}|_S \mathbf{0}]}^{\text{element A}}^T \sigma_{\Omega_k} S_{\Omega_k} + \overbrace{[\mathbf{0} \ \frac{1}{2}\mathbf{B}_{\Omega_k}|_S \mathbf{B}_{\Omega_k}|_B]}^{\text{element B}}^T \sigma_{\Omega_k} S_{\Omega_k} \quad (51) $$

or in a more effective way, as already done for stiffness matrix and strain tensor calculations:

$$ \mathbf{f}_{\text{int}} = \overbrace{[\mathbf{B}_{\Omega_k}|_A \mathbf{B}_{\Omega_k}|_S \mathbf{0}]}^{\text{element A}}^T \sigma_{\Omega_k} S_{\Omega_k} + \overbrace{[\mathbf{0} \ \mathbf{0} \ \mathbf{B}_{\Omega_k}|_B]}^{\text{element B}}^T \sigma_{\Omega_k} S_{\Omega_k}. \quad (52) $$

Equations (51) and (52) clearly show that the internal force vector can be calculated by the two elements independently since the strain-displacement matrices $\mathbf{B}_{\Omega_k}|_A$, $\mathbf{B}_{\Omega_k}|_S$, and $\mathbf{B}_{\Omega_k}|_B$ have already been evaluated during stiffness matrix evaluation (Section 5.2) and no further calculation and data exchange between the two elements are needed.

The two terms of the sum in Equation (51) can be restrained to the degrees of freedom of each element:

$$ \begin{aligned} \mathbf{f}_{\text{int}} &= \mathbf{f}_{\text{int}|A} + \mathbf{f}_{\text{int}|B}, \\ \mathbf{f}_{\text{int}|I} &= [\mathbf{B}_{\Omega_k}|_I \frac{1}{2}\mathbf{B}_{\Omega_k}|_S]^T \sigma_{\Omega_k} S_{\Omega_k} \quad \text{for } I \in \{A, B\}, \end{aligned} \quad (53) $$

where the addition symbol “+” means “assembly.”

It is worth to point out that no auxiliary element is used for internal force vector evaluation.
## 5.6 | Application of element-based FS-FEM formulation to non-tetrahedral elements and to adjacent elements of different type

Element-based FS-FEM formulation can be straightforwardly applied to any kind of finite element since it is independent of the element and smoothing domain geometries. The fact that the smoothing domain $\Omega_k$ is split in two subdomains $\Omega_k|A$ and $\Omega_k|B$ allows to compute smoothing domain related quantities inside each element ignoring the adjacent one. Elements A and B can therefore be of different type and the element-based FS-FEM formulation can be applied to meshes mixing different types of solid elements.

This is effectively another advantage of the element-based formulation since in the literature, to the best of the authors' knowledge, FS-FEM has been applied to tetrahedral meshes only.

In order to show the effectiveness of element-based FS-FEM formulation in managing a smoothing domain $\Omega_k$ extending over two elements of different type, the calculation of the stiffness matrix $K_{\Omega_k}$ for the case shown in Figure 3 will be detailed in the following section.

### 5.6.1 | Example of calculation of the stiffness matrix of a smoothing domain extending over a hexahedral and a pyramidal element

The smoothing domain $\Omega_k$ is built on the quadrangular face 5-6-7-8 in Figure 3 shared by an hexahedral and a pyramidal element. It is split in two pyramidal subdomains, one for each of the two elements sharing the quadrangular face. An auxiliary element is built by using the nodes belonging to one element only (element defined by nodes 1-2-3-4-9 in Figure 3). Matrix $K_{\Omega_k}$ is given by the sum of the contributions of hexahedral, pyramidal, and auxiliary elements (Equation 41):

$$K_{\Omega_k} = K_{\Omega_k \text{HEXA}} + K_{\Omega_k \text{PYRA}} + K_{\Omega_k \text{AUX}}, \quad (54)$$

<figure><img src="image_5.png" /><figcaption>FIGURE 3 Example of smoothing domain extending over two elements of different type. The smoothing domain is built on the quadrangular face defined by nodes 5-6-7-8. The face is shared between the hexahedral and pyramidal elements. The smoothing domain is composed by the two pyramids shaded in the figure, one inside the hexahedral element and one inside the pyramidal element, obtained by connecting the centroid of each element with the nodes defining the shared face. The auxiliary element is defined by the nodes of the two elements not belonging to the shared face.</figcaption></figure>
where the three contributions are given by:

$$ \mathbf{K}_{\Omega_k \text{HEXA}} = \begin{bmatrix} \mathbf{K}_{11}^{\Omega_k} & \cdots & \mathbf{K}_{14}^{\Omega_k} & \mathbf{K}_{15}^{\Omega_k} & \cdots & \mathbf{K}_{18}^{\Omega_k} & \mathbf{0} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{K}_{41}^{\Omega_k} & \cdots & \mathbf{K}_{44}^{\Omega_k} & \mathbf{K}_{45}^{\Omega_k} & \cdots & \mathbf{K}_{48}^{\Omega_k} & \mathbf{0} \\ \mathbf{K}_{51}^{\Omega_k} & \cdots & \mathbf{K}_{54}^{\Omega_k} & \mathbf{K}_{55}^{\Omega_k} & \cdots & \mathbf{K}_{58}^{\Omega_k} & \mathbf{0} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \ddots \\ \mathbf{K}_{81}^{\Omega_k} & \cdots & \mathbf{K}_{84}^{\Omega_k} & \mathbf{K}_{85}^{\Omega_k} & \cdots & \mathbf{K}_{88}^{\Omega_k} & \mathbf{0} \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \end{bmatrix}, \qquad (55) $$

$$ \mathbf{K}_{\Omega_k \text{PYRA}} = \begin{bmatrix} \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{K}_{59}^{\Omega_k} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{K}_{89}^{\Omega_k} \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{K}_{95}^{\Omega_k} & \cdots & \mathbf{K}_{98}^{\Omega_k} & \mathbf{K}_{99}^{\Omega_k} \end{bmatrix}, \qquad (56) $$

$$ \mathbf{K}_{\Omega_k \text{AUX}} = \begin{bmatrix} \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{K}_{19}^{\Omega_k} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{K}_{49}^{\Omega_k} \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \mathbf{K}_{91}^{\Omega_k} & \cdots & \mathbf{K}_{94}^{\Omega_k} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \end{bmatrix}. \qquad (57) $$

Different colors have been used in Equations (55)–(57) to highlight the terms belonging to each submatrix in Equation (41).

Matrices $\mathbf{K}_{il}^{\Omega_k}$ ($i, l \in [1, 9]$) are evaluated by means of Equations (36) and (39):

$$ \mathbf{K}_{il}^{\Omega_k} = (\mathbf{B}_i^{\Omega_k})^T \mathbf{C} \mathbf{B}_l^{\Omega_k}, \quad \text{for } i, l \in [1, 9]. \qquad (58) $$

The matrices $\mathbf{B}_i^{\Omega_k}$ are evaluated by means of Equations (28), (29), and (31) as follows.

First, we consider the hexahedral element. The terms $b_{ij}^{\Omega_k}$ ($i \in [1, 4]$) of matrix $\mathbf{B}_i^{\Omega_k}$ are calculated for nodes 1, 2, 3, and 4 by means of Equation (31). The calculation is done over the four triangular faces defining the boundary of smoothing domain $\Omega_k$ inside the hexahedral element only (triangles 5-6-A, 6-7-A, 7-8-A, and 8-5-A in Figure 3).

For node 9 of the pyramidal element, the same procedure applies for the evaluation of $\mathbf{B}_9^{\Omega_k}$, except that the calculations are done over the four triangular faces defining the boundary of smoothing domain $\Omega_k$ inside the pyramidal element only (triangles 5-6-B, 6-7-B, 7-8-B, and 8-5-B in Figure 3).

For shared nodes 5, 6, 7, and 8, the matrices $\mathbf{B}_i^{\Omega_k}$ ($i \in [5, 8]$) are evaluated by means of Equations (42) and (43) by adding the contributions $\mathbf{B}_i^{\Omega_k}|_{\text{HEXA}}$ and $\mathbf{B}_i^{\Omega_k}|_{\text{PYRA}}$ coming from the two adjacent elements. The calculations of these two contributions are done in exactly the same way as in the case of the calculations of $\mathbf{B}_i^{\Omega_k}$ for the nodes $i \in [1, 2, 3, 4, 9]$.

For the auxiliary element, the values of the terms $b_{ij}^{\Omega_k}$ used to calculate $\mathbf{B}_i^{\Omega_k}$ ($i \in [1, 2, 3, 4, 9]$) have already been calculated when evaluating $\mathbf{K}_{\Omega_k \text{HEXA}}$ and $\mathbf{K}_{\Omega_k \text{PYRA}}$. These values can be promptly reused to build the matrices $\mathbf{K}_{il}^{\Omega_k}$ ($i, l \in [1, 2, 3, 4, 9]$) (Equation 58) and finally $\mathbf{K}_{\Omega_k \text{AUX}}$ (Equation 57). It is worth noting that auxiliary element geometry is ignored in the calculation of $\mathbf{K}_{\Omega_k \text{AUX}}$.
It is clear that the above procedure does not take into account information concerning the element type and the same procedure applies for the two elements, hexahedral and pyramidal. The smoothing domain boundary is always decomposed in triangles, independently of the geometry of the shared face and of the two elements, and the calculations are carried out on these triangles only. Consequently, the above calculations can be carried out by the same routine in the code, independently of the type of the elements containing the smoothing domain.

## 5.7 | Shape functions for S-FEM elements

Contrary to FE method, S-FEM method does not use the derivatives of the element shape functions in its formulation ([2, 10]), as clearly shown by Equations (24)–(34). Moreover shape function values are needed only at some particular points of the boundary of the smoothing domain (Equation 31) and no explicit analytical form is required. This gives an enormous freedom in shape function definition.

Consequently, several different methods can be used to construct the shape functions provided that the following properties are satisfied at the discrete points of the element:²,⁶⁻⁸

1. Delta function: $N_i(\mathbf{x}_j) = \delta_{ij}$,

2. Partition of unity: $\sum_{i=1}^{n} N_i(\mathbf{x}) = 1$,

3. Linear compatibility: linear shape functions along element boundaries,

4. Linear consistency: $\sum_{i=1}^{n} N_i(\mathbf{x})\mathbf{x}_i = \mathbf{x}$,

5. $N_i(\mathbf{x}) \ge 0$,

where *n* is the number of nodes of the element.

In References 6-8, ES-FEM shape functions are successfully constructed by the simple averaging point interpolation method. The same method is adopted in this article for ES-FEM for the 2D cases. For FS-FEM, this method can be straightforwardly applied to tetrahedral elements because smoothing domain boundaries are composed by triangular faces only. However, for non-simplex elements, further work is required, as explained in the following section.

### 5.7.1 | Construction of FS-FEM shape functions for non-simplex elements

FS-FEM smoothing domains are built on element faces connecting face vertices with the element centroid (Figure 2). For a triangular face, the smoothing domain is given by the union of two tetrahedrons (tetrahedrons 1-2-3-G and 1-2-3-G¹ in Figure 4). For a quadrangular face, the smoothing domain is given by the union of two pyramids (pyramids 1-2-3-4-G and 1-2-3-4-G¹ in Figure 5). In the element-based formulation, the smoothing domain is processed element by element by considering only the part of the domain contained in one element. Therefore for a triangular face the two tetrahedrons composing the smoothing domain (tetrahedrons 1-2-3-G and 1-2-3-G¹ in Figure 4) are processed separately. Similarly, for a quadrangular face processing is carried out on the two pyramids separately (pyramids 1-2-3-4-G and 1-2-3-4-G¹ in Figure 5). In the following part of the section, the face generating the smoothing domain will be referred as the base face of the smoothing domain, while the other faces of the smoothing domain will be referred as the lateral faces of the domain.

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** Example of smoothing domain generated by a triangular face (face 1-2-3). Points G and G¹ are the centroids of the two elements sharing face 1-2-3. The smoothing domain is composed by the two tetrahedra 1-2-3-G and 1-2-3-G¹.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 5 Example of smoothing domain generated by a quadrangular face (face 1-2-3-4). Points G and G' are the centroids of the two elements sharing face 1-2-3. The smoothing domain is composed by the two pyramids 1-2-3-4-G and 1-2-3-4-G'.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6 Value of the shape functions at the nodes of a triangular lateral face of a smoothing domain (example of face 1-3-G in Figure 4). The numbers in parenthesis near the node/point label are the value of the shape function at the node/point. G<sub>p</sub> is the centroid of the face. (A) Example of the values when the shape function is associated to a node of the face (N<sub>1</sub>(x) in the figure). (B) Example of the values when the shape function is not associated to a node of the face.</figcaption></figure>

In order to apply the simple averaging point interpolation method, the shape function value must be known at each vertex of the smoothing domain. The values at the vertices of the base face are known since each vertex is coincident with a mesh node and delta function property 1 above applies. Shape function at element centroid G (or G<sub>1</sub>) (Figures 4 and 5), can be evaluated as follows:<sup>6-8</sup>

$$N_i(\mathbf{x}_G) = \frac{1}{n}. \qquad (59)$$

Here n is the number of nodes of the element.

In order to satisfy the linear compatibility property (property 3), linear variation of the shape function is assumed on smoothing domain faces where integration is done. Integration requires one Gauss point only and shape functions must be evaluated at the centroid of the smoothing domain face only (Equation 31). Their values can be calculated by the simple averaging point interpolation method. For a lateral face of the smoothing domain, shape function N<sub>i</sub>(x<sub>Gp</sub>) at face centroid x<sub>Gp</sub> is given by (Figure 6):

$$N_i(\mathbf{x}_{G_p}) = \frac{1 + \frac{1}{n} + 0}{3} = \frac{n+1}{3n}, \quad \text{if node } i \text{ is a vertex of the face,} \qquad (60a)$$

$$N_i(\mathbf{x}_{G_p}) = \frac{\frac{1}{n} + 0 + 0}{3} = \frac{1}{3n}, \quad \text{if node } i \text{ is not a vertex of the face.} \qquad (60b)$$

Similarly, for a triangular base face (Figure 7):

$$N_i(\mathbf{x}_{G_p}) = \frac{1+0+0}{3} = \frac{1}{3}, \quad \text{if node } i \text{ is a vertex of the face,} \qquad (61a)$$

$$N_i(\mathbf{x}_{G_p}) = \frac{0+0+0}{3} = 0, \quad \text{if node } i \text{ is not a vertex of the face.} \qquad (61b)$$

For a quadrangular base face, piecewise linear shape functions are constructed by decomposing the face into four triangles. Each triangle is obtained by connecting the face centroid G<sub>f</sub> to the vertices of each face edge, as shown in Figure 8. Integration of Equation (30) is done on each triangle by means of Equation (31). Shape function value at the
<figure><img src="image_9.png" /><figcaption>FIGURE 7 Value of the shape functions at the nodes of a triangular base face of a smoothing domain (example of face 1-2-3 in Figure 4). The numbers in parenthesis near the node label are the value of the shape function at the node. $G_p$ is the centroid of the face. (A) Example of the values when the shape function is associated to a node of the face ($N_3(\mathbf{x})$ in the figure). (B) Example of the values when the shape function is not associated to a node of the face.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8 Subdivision of a quadrangular base face of a smoothing domain in four triangles. $G_f$ is the centroid of the face.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 9 Value of the shape functions at the one of the four triangles used to split the quadrangular face in Figure 8. The numbers in parenthesis near the node label are the value of the shape function at the node/point. $G_f$ is the centroid of the quadrangular face. (A) Example of the values when the shape function is associated to a node of the face ($N_2(\mathbf{x})$ in the figure). (B) Example of the values when the shape function is not associated to a node of the face.</figcaption></figure>

quadrangle face centroid $G_f$ is calculated as in the case of element centroid G (Equation 59):

$$N_i(\mathbf{x}_{G_f}) = \frac{1}{4} \quad \text{if node } i \text{ is a vertex of the face,} \qquad (62a)$$

$$N_i(\mathbf{x}_{G_f}) = 0 \quad \text{if node } i \text{ is not a vertex of the face.} \qquad (62b)$$

In Equations (62a) and (62b), it is assumed that the shape function value of the nodes not coincident with one of the vertices of the face is zero on the face. Therefore only the four nodes defining the face contributes to the value of the shape function at $G_f$.

Simple averaging point interpolation method is then applied to each of the four triangles constituting the base face in order to calculate the shape function value at each triangle centroid $G_t$ (Figure 9) used in Equation (31):

$$N_i(\mathbf{x}_{G_t}) = \frac{0 + 1 + \frac{1}{4}}{3} = \frac{5}{12} \quad \text{if node } i \text{ is a vertex of the triangle,} \qquad (63a)$$

$$N_i(\mathbf{x}_{G_t}) = \frac{0 + 0 + \frac{1}{4}}{3} = \frac{1}{12} \quad \text{if node } i \text{ is not a vertex of the triangle.} \qquad (63b)$$
5.7.2 | Standard patch test

The solid patch test proposed by Macneal and Harder¹⁵ is used to test the non-simplex formulation of FS-FEM. Further analysis of the sensitivity to mesh distortion will be done later in Section 6.

In the solid patch test, the following prescribed boundary conditions are applied on the external faces of a unit cube (Figure 10):

$$ 
\begin{aligned}
u(\mathbf{x}) &= 10^{-3} (2x + y + z) / 2, \\
v(\mathbf{x}) &= 10^{-3} (x + 2y + z) / 2, \\
w(\mathbf{x}) &= 10^{-3} (x + y + 2z) / 2,
\end{aligned} 
\qquad (64) $$

where $u(\mathbf{x})$, $v(\mathbf{x})$, and $w(\mathbf{x})$ are the displacements in x, y, and z directions, respectively, at a point of coordinates $\mathbf{x} = [x, y, z]^T$.

The length of the edges of the unit cube is taken equal to 1 m. The cube material is assumed linear elastic and isotropic. Its Young's modulus and Poisson's ratio are $E = 1$ MPa and $\nu = 0.25$, respectively.

The theoretical solution for stress and strain fields is given by:¹⁵

$$ \epsilon_x = \epsilon_y = \epsilon_z = \gamma_{xy} = \gamma_{yz} = \gamma_{zx} = 10^{-3}, \qquad (65a) $$

$$ 
\begin{aligned}
\sigma_x &= \sigma_y = \sigma_z = 2000 \text{ Pa}, \\
\tau_{xy} &= \tau_{yz} = \tau_{zx} = 400 \text{ Pa}.
\end{aligned} 
\qquad (65b) $$

The hexahedral mesh used in the solid patch test is shown in Figure 10. The coordinates of the internal nodes are given in Table 1.

The comparison between the theoretical solution and the FS-FEM solution is done at the internal nodes given in Table 1 for the displacements and at each smoothing cell for the stress and strain fields. This choice is based on the fact that the stress and strain fields given by the theoretical solution are constant in space (Equations 65a and 65b) and the stress and strain fields calculated by FS-FEM are constant over a smoothing domain.

The comparison between the theoretical and numerical values of the displacement is done by evaluating the relative error for each displacement component:

$$ 
\begin{aligned}
e_u(\mathbf{x}_i) &= \frac{|u(\mathbf{x}_i) - u_{\text{SFEM}}(\mathbf{x}_i)|}{u(\mathbf{x}_i)}, \\
e_v(\mathbf{x}_i) &= \frac{|v(\mathbf{x}_i) - v_{\text{SFEM}}(\mathbf{x}_i)|}{v(\mathbf{x}_i)}, \\
e_w(\mathbf{x}_i) &= \frac{|w(\mathbf{x}_i) - w_{\text{SFEM}}(\mathbf{x}_i)|}{w(\mathbf{x}_i)},
\end{aligned} 
\qquad (66) $$

<figure><img src="image_12.png" /><figcaption>FIGURE 10 Unit cube patch test. Internal node number is shown.</figcaption></figure>
<table><caption>TABLE 1 Coordinates of the internal nodes of the unit cube patch test¹⁵</caption><thead><tr><th>Node number</th><th>x</th><th>y</th><th>z</th></tr></thead><tbody><tr><td>1</td><td>0.249</td><td>0.342</td><td>0.192</td></tr><tr><td>2</td><td>0.826</td><td>0.288</td><td>0.288</td></tr><tr><td>3</td><td>0.850</td><td>0.649</td><td>0.263</td></tr><tr><td>4</td><td>0.273</td><td>0.750</td><td>0.230</td></tr><tr><td>5</td><td>0.320</td><td>0.186</td><td>0.643</td></tr><tr><td>6</td><td>0.677</td><td>0.305</td><td>0.683</td></tr><tr><td>7</td><td>0.788</td><td>0.693</td><td>0.644</td></tr><tr><td>8</td><td>0.165</td><td>0.745</td><td>0.702</td></tr></tbody></table>

<table><caption>TABLE 2 Maximum value of the errors for the unit cube patch test</caption><thead><tr><th></th><th>e<sub>u</sub>(x<sub>i</sub>)</th><th>e<sub>v</sub>(x<sub>i</sub>)</th><th>e<sub>w</sub>(x<sub>i</sub>)</th><th>e<sub>σ</sub><sup>Ω<sub>k</sub></sup></th><th>e<sub>ε</sub><sup>Ω<sub>k</sub></sup></th></tr></thead><tbody><tr><td>Maximum value</td><td>0.0</td><td>0.0</td><td>0.0</td><td>2.56e-15</td><td>1.08e-15</td></tr></tbody></table>

Note: Errors on displacements components ($e_u(x_i)$, $e_v(x_i)$, $e_w(x_i)$) are evaluated at mesh nodes. Errors on stress and strain tensor components ($e_σ^{Ω_k}$, $e_ε^{Ω_k}$) are evaluated at each smoothing domain.

where $e_u(x_i)$, $e_v(x_i)$, and $e_w(x_i)$ are the relative errors at node i (with coordinates $x_i$ given in Table 1) on x, y, and z displacements, respectively, $u(x_i)$, $v(x_i)$, and $w(x_i)$ are the theoretical displacements at node i given by Equation (64) and $u_{SFEM}(x_i)$, $v_{SFEM}(x_i)$, and $w_{SFEM}(x_i)$ are the corresponding S-FEM numerical values.

The error for the stress and strain fields is given by the maximum value of the relative error calculated on each component of the stress and strain tensors at each smoothing domain $Ω_k$. For the stress field, the error reads:

$$e_{\sigma}^{\Omega_k} = \max_{\sigma_{\text{cmp}} \in A} \left( \frac{|\sigma_{\text{cmp}} - \sigma_{\text{cmp}}^{\Omega_k}|}{\sigma_{\text{cmp}}} \right), \quad \text{for } A = \{\sigma_x, \sigma_y, \sigma_z, \tau_{xy}, \tau_{yz}, \tau_{xz}\}, \qquad (67)$$

where $e_σ^{Ω_k}$ is the maximum relative error on the stress tensor components at smoothing domain $Ω_k$, $σ_{cmp}$ is one component of the stress tensor of the theoretical solution given by Equation (65b) and $σ_{cmp}^{Ω_k}$ is the corresponding S-FEM numerical value.

For the strain field, the error reads:

$$e_{\epsilon}^{\Omega_k} = \max_{\epsilon_{\text{cmp}} \in A} \left( \frac{|\epsilon_{\text{cmp}} - \epsilon_{\text{cmp}}^{\Omega_k}|}{\epsilon_{\text{cmp}}} \right), \quad \text{for } A = \{\epsilon_x, \epsilon_y, \epsilon_z, \gamma_{xy}, \gamma_{yz}, \gamma_{xz}\}, \qquad (68)$$

where $e_σ^{Ω_k}$ is the maximum relative error on the stress tensor components at smoothing domain $Ω_k$, $ε_{cmp}$ is one component of the stress tensor of the theoretical solution given by Equation (65a) and $ε_{cmp}^{Ω_k}$ is the corresponding S-FEM numerical value.

The errors are given in Table 2. They are of the same order of magnitude of the machine precision. It can be concluded that the FS-FEM formulation used for non-simplex elements passes the patch test and is able to represent a constant stress and strain state.

# 6 | NUMERICAL EXPERIMENTS

In order to verify that the proposed element-based formulation for ES-FEM and FS-FEM methods works, numerical experiments have been carried out for 2D and 3D cases. All the elements commonly used in finite element analysis have
been tested and mesh convergence has been verified by comparing numerical solution obtained by the proposed approach with finite element and reference solutions available in the literature.

Convergence tests have been carried out first on regular meshes with element aspect ratios equal or close to 1. As done in Reference 10, in order to check the effect of element distortion on the stability and accuracy of the proposed formulation, further convergence tests have been done on distorted meshes obtained from the regular one by moving the nodes as follows:

$$ 
\begin{aligned}
x &= x_0 + \alpha_{ir} \cdot r_x \cdot \Delta h_x, \\
y &= y_0 + \alpha_{ir} \cdot r_y \cdot \Delta h_y, \\
z &= z_0 + \alpha_{ir} \cdot r_z \cdot \Delta h_z,
\end{aligned} 
\qquad (69) $$

where $x_0 = [x_0, y_0, z_0]^T$ and $x = [x, y, z]^T$ are the coordinates of the node before and after the distortion, respectively, $\alpha_{ir} \in [0, 1]$ is an irregularity factor (element is not distorted for $\alpha_{ir} = 0$ and distortion increases with $\alpha_{ir}$), $r_x, r_y$, and $r_z$ are three computer generated random numbers included in the interval $[-0.5, 0.5]$ and $\Delta h_x, \Delta h_y$, and $\Delta h_z$ are the mean edge size of the finite elements in the mesh in $x, y$, and $z$ directions, respectively.

The element-based formulation for ES-FEM and FS-FEM methods have been implemented in Code_Aster, a free finite element software developed by EDF R&D for solid mechanics.16 All the numerical results shown in this article have been obtained by means of this modified version of Code_Aster. Element matrix computation routines of existing finite elements have been replaced by a generic S-FEM routine. Auxiliary elements are created dynamically during the simulation as “late elements” (“éléments tardifs” in the original French documentation of the code). These elements are standard elements available in the code library. Any routine in the code can dynamically add them to an existing mesh provided by the user. They do not belong to the mesh and they are not shown in calculation results, but they participate in code calculations, like stiffness matrices calculations.

## 6.1 | 2D test case

For the 2D case, the elastic case of the cantilever beam shown in Figure 11 has been studied. The thickness of the beam is assumed to be extremely small with respect to its length $L$ and its height $D$. The beam is therefore in a plane stress condition. The parabolic distributed shear force applied on the free edge on the right (see Figure 11) is given by:

$$ t_y(y) = \frac{6P}{D^3}y^2 - \frac{3P}{2D}, \qquad (70) $$

where $D$ is the height of the beam (see Figure 11) and $P$ is the resultant shear force applied on the free edge.

The material is assumed linear elastic and isotropic. Its Young’s modulus and Poisson’s ratio are $E$ and $\nu$, respectively.

The analytical solution of the problem is available in Reference 6. The displacement field is given by the following equations:

$$ 
\begin{aligned}
u_x &= \frac{Py}{6EI} \left[ (6L - 3x)x - (2 + \nu) \left(y^2 - \frac{D^2}{4}\right) \right], \\
u_y &= -\frac{P}{6EI} \left[ 3\nu y^2(L - x) + (4 + 5\nu) \frac{D^2x}{4} + (3L - x)x^2 \right],
\end{aligned} 
\qquad (71) $$

where $I = \frac{D^3}{12}$ is the area moment of inertia per unit width of the beam cross-section.

<figure><img src="image_13.png" /><figcaption>FIGURE 11 Cantilever beam subjected to a parabolic shear force on the free edge on the right</figcaption></figure>
The plane stress field is given by the following equations:

$$ 
\begin{aligned}
\sigma_x &= \frac{P(L-x)y}{I}, \\
\sigma_y &= 0, \\
\tau_{xy} &= -\frac{P}{2I} \left( \frac{D^2}{4} - y^2 \right).
\end{aligned} 
\qquad (72) $$

For the numerical test, the problem parameters have been taken as $E = 2.1 \times 10^{11}$ Pa, $\nu = 0.3$, $D = 50$ m, $L = 500$ m, and $P = 10$ kN.⁶ The thickness of the beam has been fixed to 1 m. First, two uniform meshes of four-node quadrilateral elements (Q4) and linear triangular elements (T3) have been used. The T3 mesh has been obtained from the Q4 mesh by splitting each quadrilateral element in two triangles using one of its diagonals (Figure 12). The aspect ratio of the Q4 elements has been fixed to 1 to ensure high precision of the finite element solution. Distorted meshes have been generated from these regular meshes by moving the mesh nodes by means of Equation (69). An examples of distorted Q4 mesh is shown in Figure 13 ($\alpha_{ir} = 0.7$). The highest irregularity factor used in the tests is $\alpha_{ir} = 0.7$, which corresponds to a very distorted mesh. For higher irregularity factors, overlapping elements appear in the mesh (Figure 14), invalidating the hypothesis at the base of S-FEM theory (nonoverlapping smoothing domains, see section 4).

For the regular mesh case ($\alpha_{ir} = 0$), mesh convergence is shown in Figure 15 for the FEM solution and for the element-based S-FEM formulation presented in this article. Solutions are compared in terms of elastic strain energy. In Figure 16, the comparison is done in terms of the relative error in energy norm, calculated as follows:

$$ \text{Error} = \sqrt{\frac{\int (\boldsymbol{\sigma}_{\text{FE}} - \boldsymbol{\sigma}_{\text{ref}}) : \boldsymbol{C}^{-1} : (\boldsymbol{\sigma}_{\text{FE}} - \boldsymbol{\sigma}_{\text{ref}})}{\int \boldsymbol{\sigma}_{\text{ref}} : \boldsymbol{C}^{-1} : \boldsymbol{\sigma}_{\text{ref}}}}, \qquad (73) $$

where $C$ is the stiffness tensor of the generalized Hooke's law (Equation 9) and $\sigma_{FE}$ and $\sigma_{ref}$ are the finite element (FEM or S-FEM) and reference (Equation 72) stress tensors, respectively. Integration in Equation (73) is carried out over the whole mesh.

The comparison between FEM and ES-FEM solutions for the distorted mesh case ($\alpha_{ir} = 0.7$) is shown in Figures 17 and 18 for the elastic strain energy and relative error in energy norm, respectively.

<figure><img src="image_14.png" /><figcaption>FIGURE 12 Examples of Q4 and T3 meshes used for the cantilever beam of Figure 11. T3 mesh is obtained from Q4 mesh by splitting each quadrilateral element in two triangular elements.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 13 Example of a distorted Q4 mesh ($\alpha_{ir} = 0.7$)</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIGURE 14 Example of overlapping T3 elements generated by high irregularity factors</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 15 Results for test case of Figure 11. Comparison between the elastic strain energy obtained from S-FEM and FEM models for different element sizes and the value from the analytical solution for T3 and Q4 meshes ($\alpha_{ir} = 0$)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 16 Results for test case of Figure 11. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for T3 and Q4 meshes ($\alpha_{ir} = 0$)</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIGURE 17 Results for test case of Figure 11. Comparison between the elastic strain energy obtained from S-FEM and FEM models for different element sizes and the value from the analytical solution for T3 and Q4 meshes ($\alpha_{ir}$ = 0.7)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 18 Results for test case of Figure 11. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for T3 and Q4 meshes ($\alpha_{ir}$ = 0.7)</figcaption></figure>

As expected,⁶ the element-based ES-FEM method converges toward the reference solution with mesh refinement and the accuracy and mesh convergence rate are always higher than for the FEM, independently of the distortion of the mesh. These results show that ES-FEM is less sensitive to mesh distortion than FEM.

## 6.2 | 3D test case

The cubic cantilever shown in Figure 19 is subjected to a uniform pressure of 1 Pa on the upper face.¹⁰ The face lying in the y-z plane ($x = 0$) is fixed. The material is linear elastic (Young's modulus $E = 1$ Pa, Poisson's ratio $\nu = 0.25$). A reference solution is given in Reference 10. It has been calculated by means of a FE model with a very refined mesh of 10-node tetrahedral elements. The reference displacement at point A at (1.0, 1.0, -0.5) (Figure 19) is 3.3912 m, the reference strain energy is 0.9486 N m. The theoretical solution is not available. Therefore, in order to calculate the relative error in energy norm of the element-based FS-FEM solution by means of Equation (73), a reference stress field is calculated by means of a FE model using a refined regular mesh. The mesh is composed by 64,000 cubic 8-node hexahedral elements (74,080 nodes). The displacement at point A obtained by this refined FE model is 3.3910 m, the strain energy is 0.948633 N m. These results are in excellent agreement with the solution given in Reference 10. Relative error in strain energy norm (Equation 73) is evaluated on this refined mesh using the calculated reference stress field and the element-based FS-FEM stress field from the coarser mesh projected on this refined mesh.
<figure><img src="image_21.png" /><figcaption>**FIGURE 19** Cubic cantilever subjected to a uniform pressure of 1 Pa on the upper free face</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**FIGURE 20** Example of uniform hexahedral mesh of cubic cantilever of Figure 19</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**FIGURE 21** Example of the mesh of uniform wedge elements obtained by splitting the hexahedral mesh of Figure 20. Wireframe view of hexahedral mesh of Figure 20 is shown. Wedge elements are shrunk to show how they have been obtained from hexahedral elements.</figcaption></figure>

Several simulations have been carried out in order to test mesh convergence of the element-based FS-FEM formulation and its applicability to any kind of standard element used in FE simulations. First, regular meshes have been tested. The eight-node hexahedral element (H8) mesh is shown in Figure 20. The element aspect ratio is 1. The mesh is refined by regular element subdivision preserving the element aspect ratio. Meshes using different element types have been built from the hexahedral element mesh by element subdivision: six-node wedge element (W6) mesh, four-node tetrahedral element (T4) mesh and five-node pyramidal element (P5) mesh are shown in Figures 21–23, respectively. Distorted meshes have been generated from these regular meshes by moving the mesh nodes by means of Equation (69). An example of distorted H8 mesh is shown in Figure 24. The highest irregularity factor used in the tests is $\alpha_{ir} = 0.7$, which corresponds to a very distorted mesh. For higher irregularity factors, the FE software is not able to calculate equivalent nodal forces for the distributed pressure load applied to the top surface of the cube due to the high distortion of the surface mesh. Moreover, for higher irregularity factors, some overlapping elements appear in the mesh, as in the 2D test case described in Section 6.1.
<figure><img src="image_24.png" /><figcaption>**FIGURE 22** Example of the mesh of uniform tetrahedral elements obtained by splitting the hexahedral mesh of Figure 20. Wireframe view of hexahedral mesh of Figure 20 is shown. Tetrahedral elements are shrunk to show how they have been obtained from hexahedral elements.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**FIGURE 23** Example of the mesh of uniform pyramidal elements obtained by splitting the hexahedral mesh of Figure 20. Wireframe view of hexahedral mesh of Figure 20 is shown. Pyramidal elements are shrunk to show how they have been obtained from hexahedral elements.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**FIGURE 24** Example of a distorted H8 mesh ($\alpha_{ir}$ = 0.7)</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>FIGURE 25 Results for test case of Figure 19. Comparison between the values of the displacement of point A in Figure 19 obtained from the FS-FEM and FEM models for meshes of different element type and different element size ($\alpha_{ir} = 0$). A reference solution from Reference 10 is shown as well.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIGURE 26 Results for test case of Figure 19. Comparison between the values of the strain energy obtained from the FS-FEM and FEM models for meshes of different element type and different element size ($\alpha_{ir} = 0$). A reference solution from Reference 10 is shown as well.</figcaption></figure>

The comparison between reference, finite element and element-based FS-FEM solutions for the regular meshes is shown in Figures 25 and 26 for displacement of point A and strain energy, respectively. The relative error in energy norm is shown in Figures 27–30. As expected,¹⁰ the element-based FS-FEM method converges toward the reference solution with mesh refinement and the mesh convergence rate is higher than for the FEM. The accuracy is always higher than for the FEM.

A similar comparison is done for the distorted meshes (irregularity factor $\alpha_{ir} = 0.7$) in Figures 31–36. As expected, the error is greater for the distorted meshes than for the regular meshes. The S-FEM accuracy and convergence rate with mesh refinement are always higher than for the FEM. It should be pointed out that for the distorted meshes the FE software complains about elements in the mesh with nonconstant sign of the Jacobian of the isoparametric transformation matrix, which is not the case for S-FEM since no isoparametric transformation is used (Section 4). Variation of Jacobian sign over an element is detrimental for the numerical solution accuracy, as shown for example in Figure 37 for a hexahedral mesh. These results show that S-FEM is less sensitive to mesh distortion than FEM.

## 6.3 | Computational cost of ES-FEM and FS-FEM methods

As shown in References 2, 7, 8, and 10, stiffness matrices generated by ES-FEM and FS-FEM methods are less coarse than the ones generated by FEM. Consequently, the computational cost for stiffness matrix inversion is much higher for ES-FEM and FS-FEM than for FEM. However, solution accuracy is higher for ES-FEM and FS-FEM than for FEM, as
<figure><img src="image_29.png" /><figcaption>FIGURE 27 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for H8 meshes ($\alpha_{ir} = 0.0$)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE 28 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for W6 meshes ($\alpha_{ir} = 0.0$)</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>FIGURE 29 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for T4 meshes ($\alpha_{ir} = 0.0$)</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>FIGURE 30 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for P5 meshes ($\alpha_{ir} = 0.0$)</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>FIGURE 31 Results for test case of Figure 19. Comparison between the values of the displacement of point A in Figure 19 obtained from the FS-FEM and FEM models for meshes of different element type and different element size ($\alpha_{ir} = 0.7$). A reference solution from Reference 10 is shown as well.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>FIGURE 32 Results for test case of Figure 19. Comparison between the values of the strain energy obtained from the FS-FEM and FEM models for meshes of different element type and different element size ($\alpha_{ir} = 0.7$). A reference solution from Reference 10 is shown as well.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>FIGURE 33 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for H8 meshes ($\alpha_{ir} = 0.7$)</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>FIGURE 34 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for W6 meshes ($\alpha_{ir} = 0.7$)</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>FIGURE 35 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for T4 meshes ($\alpha_{ir} = 0.7$)</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>FIGURE 36 Results for test case of Figure 19. Relative error in energy norm as defined in Equation (73) for the S-FEM and FEM numerical solutions for different element sizes for P5 meshes ($\alpha_{ir} = 0.7$)</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>FIGURE 37 Effect of the change in the sign of the Jacobian on the FEM stress field, $\sigma_y$ stress component (Pa) for the external face at $z = 0.5$ of the cube of Figure 19. The H8 mesh is distorted with an irregularity factor $\alpha_{ir} = 0.7$. On the left: FEM solution. On the right: FS-FEM solution</figcaption></figure>

shown in Sections 6.1 and 6.2. Therefore, a more relevant parameter for the performance comparison of the methods is the computational efficiency, as done in References 2 and 7-10: computation time should be compared at a fixed accuracy rather than a fixed mesh. Since FEM accuracy is lower than ES-FEM and FS-FEM accuracy, in order to get the same accuracy, a finer mesh is needed for FEM compared to ES-FEM and FS-FEM, which increase the computation time.

The element-based formulation presented in this article adds an overhead to the standard S-FEM formulation because some time is needed to create and manage auxiliary elements. In order to evaluate this overhead, an example for the hexahedral mesh from the 3D case presented in Section 6.2 will be used. In Figure 27, the relative error in energy norm is shown for this case. For FEM, the biggest case is the one with 9281 nodes (the last point on the right on the FEM curve in the figure) for which the error is close to $9 \times 10^{-2}$. The same error is obtained by FS-FEM with a coarser mesh composed by 3375 nodes. Computation times for these two models are shown in Table 3, where the model preparation time and the solver time are shown. The model preparation time includes everything is needed to prepare the model for the computation, whereas the solver time is the time needed to solve the system of linear algebraic equations produced by the numerical discretization. For the FS-FEM method, the model preparation time includes also the time needed to create and manage the auxiliary elements necessary to the element-based formulation.

As clearly shown in Table 3, for a given mesh, solver time for FEM is lower than for FS-FEM. However, for a given accuracy, solver time for FEM (9281 nodes mesh) is higher than for FS-FEM (3375 nodes mesh). Model preparation time
<table><caption>TABLE 3 Computation time for FEM and FS-FEM for two different hexahedral meshes (element aspect ratio equal to 1)</caption><thead><tr><th rowspan="2">Mesh nodes</th><th colspan="2">FEM</th><th colspan="2">FS-FEM</th></tr><tr><th>Preparation (s)</th><th>Solver (s)</th><th>Preparation (s)</th><th>Solver (s)</th></tr></thead><tbody><tr><td>3375</td><td>0.07</td><td>1.10</td><td>0.17</td><td>2.94</td></tr><tr><td>9281</td><td>0.11</td><td>3.95</td><td>0.38</td><td>14.59</td></tr></tbody></table>

*Note*: “Preparation” means the preparation computation time needed to prepare the model for the simulation, including the time needed to create and manage auxiliary elements for FS-FEM. “Solver” means the time needed by the linear solver to solve the discretized model.

is much lower than solver time both for FEM and FS-FEM. An estimation of the overhead of element-based FS-FEM formulation can be obtained by comparing the model preparation time for FEM and FS-FEM for a given mesh: for the two meshes shown in Table 3, model preparation time for FS-FEM is about 3 times the model preparation time for FEM. However, this overhead is much smaller than the solver time. Therefore it does not impact the total computation time significantly.

## 7 | CONCLUSIONS

An element-based formulation for ES-FEM and FS-FEM models has been presented in detail. It has been shown that S-FEM smoothing domain-based equations can be equivalently rewritten on an element basis if an auxiliary element is introduced in the mesh for each smoothing domain. This auxiliary element is used for stiffness matrix calculations only. The element-based formulation allows the implementation of ES-FEM and FS-FEM models in an existing finite element code for any type of finite element because it fits into the existing element-based architecture of FE codes.

The effectiveness of the element-based formulation has been shown by numerical tests and comparison with existing solutions in literature. Finally, the FS-FEM model has been extended to non-tetrahedral elements by formulating shape functions for non-simplex smoothing domains. The ability of the element-based formulation to manage any type of solid element has been shown by numerical examples.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ORCID

Daniele Colombo https://orcid.org/0000-0001-7249-9486

Ralf Frotschner https://orcid.org/0000-0002-7313-7353

Manfred Staat https://orcid.org/0000-0003-4363-6570

## REFERENCES

1. Zeng W, Liu GR. Smoothed finite element methods (S-FEM): an overview and recent developments. *Arch Comput Methods Eng.* 2018;25:397-435. doi:10.1007/s11831-016-9202-3

2. Liu GR, Thoi TN. Smoothed Finite Element Method. CRC Press; 2010.

3. Liu GR, Nguyen-Xuan H, Nguyen-Thoi T. A theoretical study on the smoothed FEM (S-FEM) models: properties, accuracy and convergence rates. *Int J Numer Methods Eng.* 2010;84:1222-1256. doi:10.1002/nme.2941

4. Liu GR, Dai KY, Nguyen TT. A smoothed finite element method for mechanics problems. *Comput Mech.* 2007;39:859-877. doi:10.1007/s00466-006-0075-4

5. Liu GR, Nguyen-Thoi T, Nguyen-Xuan H, Lam KY. A node-based smoothed finite element method (NS-FEM) for upper bound solutions to solid mechanics problems. *Comput Struct.* 2009;87:14-26. doi:10.1016/j.compstruc.2008.09.003

6. Liu GR, Nguyen-Thoi T, Lam KY. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. *J Sound Vib.* 2009;320:1100-1130. doi:10.1016/j.jsv.2008.08.027

7. Dai KY, Liu GR, Nguyen TT. An n-sided polygonal smoothed finite element method (nSFEM) for solid mechanics. *Finite Elem Anal Des.* 2007;43:847-860. doi:10.1016/j.finel.2007.05.009

8. Nguyen-Thoi T, Liu GR, Nguyen-Xuan H. An n-sided polygonal edge-based smoothed finite element method (nES-FEM) for solid mechanics. *Int J Numer Methods Biomed Eng.* 2010;27(9):1446-1472. doi:10.1002/cnm.1375

9. Cui X, Liu GR, Li G, Zhang G, Zheng G. Analysis of plates and shells using an edge-based smoothed finite element method. *Comput Mech.* 2010;45:141-156. doi:10.1007/s00466-009-0429-9
10. Nguyen-Thoi T, Liu GR, Lam KY, Zhang GY. A face-based smoothed finite element method (FS-FEM) for 3D linear and geometrically non-linear solid mechanics problems using 4-node tetrahedral elements. *Int J Numer Methods Eng.* 2009;78:324-353. doi:10.1002/nme.2491

11. Bathe KJ. *Finite Element Procedures*. Prentice-Hall; 1996.

12. Zienkiewicz OC, Taylor RL. *The Finite Element Method*. Vol 1. 5th ed. Butterworth-Heinemann; 2000.

13. Cook RD. *Finite Element Modeling for Stress Analysis*. John Wiley & Sons, Inc.; 1995.

14. Liu GR. A G space theory and a weakened weak ($W^2$) form for a unified formulation of compatible and incompatible methods: Part II applications to solid mechanics problems. *Int J Numer Methods Eng.* 2009;81(9):1127-1156. doi:10.1002/nme.2720

15. Macneal RH, Harder RL. A proposed standard set of problems to test finite element accuracy. *Finite Elem Anal Des.* 1985;1:3-20.

16. EDF R&D, code_aster finite element code. Accessed December 23, 2021. http://www.code-aster.org

**How to cite this article:** Colombo D, Drira S, Frotscher R, Staat M. An element-based formulation for ES-FEM and FS-FEM models for implementation in standard solid mechanics finite element codes for 2D and 3D static analysis. *Int J Numer Methods Eng.* 2023;124(2):402-433. doi: 10.1002/nme.7126