Monolithic coupling of implicit material point method with finite element method

Three different domain configurations are used in the mathematical formulation. The reason why the configurations are three is made explicit in the following sections. Let $T>0$ and let $\Omega$ be a given open and bounded domain such that, for all $t\in[0,T]$, $\Omega(t)\subset\mathbb{R}^{d}$. Then, the undeformed configuration of $\Omega$ is indicated with the notation $\widehat{\Omega}$ and it represents the position of the domain in its original state at time $t=0$. For the deformed configuration $\Omega(t)$, no superscript is used, and the notation $\Omega(t)$ just refers to the domain in its deformed configuration at time $t>0$. This configuration is obtained from the undeformed configuration with a translation. Specifically, if $\widehat{\bm{x}}$ represents any point of $\widehat{\Omega}$, then

where $\bm{u}(\widehat{\bm{x}},t)$ is the displacement field. The third configuration to be introduced is the reference configuration, denoted by $\widetilde{\Omega}$. We remark that, given $n\geq 0$, such a configuration does not change for any $t\in(t^{n},t^{n+1}]$. The reference configuration can also be obtained from the undeformed configuration with a translation, as follows

where the field $\bm{u}^{reset}$ is referred to as the reset displacement field associated with the given time interval $(t^{n},t^{n+1}]$. The definition of such a field is given in the next section. Any point in the deformed configuration can be also obtained from the reference configuration as

where the field $\widetilde{\bm{u}}$ is the displacement field with respect to the reference configuration, which satisfies

Schematics for the three domain configurations can be found in Figure 1.

The general considerations about domain configurations are applied to the MPM and FEM domains as explained next.

Let $T>0$ and $[0,T]\subset\mathbb{R}$. For all $t\in[0,T]$, let $\Omega_{mpm}(t)$ and $\Omega_{b}(t)$ be open and bounded subsets of $\mathbb{R}^{d}$. Assume that $\Omega_{mpm}(t)\subset\Omega_{b}(t)$ for all $t\in[0,T]$. From now on, $\Omega_{mpm}(t)$ represents the body discretized with the MPM particles. The set $\Omega_{b}(t)$ is discretized with a regular finite element grid $\mathcal{T}_{b}$ [14, 8], used as a background grid for the MPM body. While $\Omega_{mpm}(t)$ may be subject to roto-translation and large deformations due to the movement of the MPM body, $\Omega_{b}(t)$ undergoes only limited deformations, because it follows the MPM body motion only in the time interval $(t^{n},t^{n+1}]$. At the beginning of each time step, $\Omega_{b}(t)$ is reset to its reference configuration $\widetilde{\Omega}_{b}$, as it will be explained in the next section.

Denoting by $\mathcal{E}_{p}$ the element of the background grid that hosts a given particle $p$, we introduce

where $t\in(t^{n},t^{n+1}]$, and $\mathcal{T}_{ab}\subset\mathcal{T}_{b}$ is the active background grid, composed of all elements $\mathcal{E}$ of $\mathcal{T}_{b}$ that host at least one particle $p$. Similarly, we define $\mathcal{T}_{ib}\equiv\mathcal{T}_{b}\setminus\mathcal{T}_{ab}$ to be the inactive background grid, composed of all elements of $\mathcal{T}_{b}$ that do not enclose any particle. Note that the time dependence of $\Omega_{ab}(t)$ is due to the time dependence of $\mathcal{E}_{p}$ (for given $p$) in the time interval $(t^{n},t^{n+1}]$. Moreover, because the background grid is reset at the beginning of each time step, a particle $p$ may leave a particular element $\mathcal{E}_{p}$ and move into one of its neighboring elements, so the elements of $\mathcal{T}_{ab}$ and $\mathcal{T}_{ib}$ may also change at each time step. Schematics of the grids here described are shown in Figure 2. By definition, $\Omega_{ab}(t)$ has the property that $\Omega_{mpm}(t)\subseteq\Omega_{ab}(t)\subseteq\Omega_{b}(t)$. The set $\Omega_{ab}(t)$, associated with the MPM active background grid, is what we refer to as the MPM domain. Such a choice is motivated by the fact that the MPM momentum equation is solved only on the active background grid, and its elements consequently undergo deformations that modify their reference configuration.

Concerning the FEM, since the triangulation entirely covers the body to discretize, the open and bounded set $\Omega_{fem}(t)\subset\mathbb{R}^{d}$ is used to characterize both the body discretized with FEM and the FEM domain. The interfaces between the FEM domain and the background grid $\Omega_{b}(t)$ and between the FEM domain and the active background grid $\Omega_{ab}(t)$ that arise during the coupling procedure are indicated with the letters $\Gamma$ and $\Gamma_{a}$, respectively.

In the classical MPM method, where no coupling between MPM and FEM occurs, at the beginning of each time step interval $(t^{n},t^{n+1}]$, the set $\Omega_{b}(t)$ is reset to its undeformed configuration, $\widehat{\Omega}_{b}$, to avoid mesh entanglement. Namely, for any $t\in(t^{n},t^{n+1}]$, if $I$ refers to any node of the grid, the grid displacement $\bm{u}_{I}$ associated to the node $I$ is such that

Moreover, in general

and as a result $\bm{u}_{I}(t)$ is discontinuous in time. On the other hand, in the classical uncoupled FEM method, the grid follows the solid body deformation and its displacement is continuous in time

With a monolithic coupling approach, there is a unique displacement field defined on the MPM background grid and on the finite element grid discretizing the FEM body. Hence, when the coupling between MPM and FEM is considered, the displacement at the grid points of $\mathcal{T}_{b}$ that lie on the interface $\Gamma$ between $\Omega_{b}(t)$ and the FEM body cannot be reset to zero. This implies that at every instant of time, the deformed configuration of the MPM domain is composed of elements that are deformed, not only because of the motion of the MPM body, but also because of the FEM body motion.

The requirement of resetting the MPM background grid to its undeformed configuration at the beginning of every time step is what motivates the introduction of the reset displacement field to be used for the coupled problem. More specifically, when MPM is coupled with FEM, the background grid has to fulfill the following constraints: in the FEM domain and on the interface $\Gamma$, it has to follow the solid deformation; in the MPM domain, in order to enforce domain continuity, it has to follow the deformation of $\Gamma$, coming from the FEM domain. Moreover, to avoid mesh entanglement, it is desirable to reset the background grid to a configuration similar to the undeformed configuration on all MPM nodes that are sufficiently far from $\Gamma$.

Because in the applications described in this paper the FEM solid body is subject to small deformations, the background grid constraints are satisfied with the following definition of the reset displacement field $\bm{u}_{I}^{reset}$. Let $\mathcal{I}$ be the set of grid points of $\mathcal{T}=\mathcal{T}_{b}\cup\mathcal{T}_{fem}$, where $\mathcal{T}_{fem}$ denotes the finite element grid used to discretize the FEM body. For all $t\in(t^{n},t^{n+1}]$, the reset displacement field satisfies

where $\bm{x}_{I}$ denotes the position of node $I\in\mathcal{I}$. When MPM is coupled with FEM, Eq. (5) is replaced with

By making this choice, continuity of the domain is preserved across $\Gamma$, and the grid $\mathcal{T}_{b}$ follows the solid body deformation on $\Omega_{fem}(t)$, resulting in a continuous displacement field $\bm{u}_{I}(t)$. As a consequence of Eq. (9), at the beginning of each time step, the deformed configuration of the background grid coincides with its reference configuration, although it is in general different from its undeformed configuration $\widehat{\Omega}_{b}$. This is more perceptible on those elements close to $\Gamma$ that may deform more to follow the FEM body motion. If the FEM body deformation is small, than the deformation of the background grid will also be small. To fulfill the classical MPM requirement in Eq. (5), a new MPM displacement field based on the reference background grid configuration is defined as

By Eq. (9), this new field satisfies

as in the classical uncoupled MPM method. The reference configuration is the configuration of the MPM domain associated with the field $\widetilde{\bm{u}}_{I}(t)$.

For ease of notation, during the rest of this paper, the time dependence of $\Omega_{mpm}$, $\Omega_{ab}$, and of the fields involved in the formulation is not made explicit. The mass conservation is automatically satisfied by the material point method [41], therefore the governing equations for the MPM body consist of momentum conservation with appropriate boundary and initial conditions. For simplicity, we consider zero boundary conditions for displacement and normal stress

The case of Robin and mixed boundary conditions is not addressed at present. In Eq.(11), $\bm{u}$ represents the displacement, $\bm{\sigma}$ is the Cauchy stress tensor, $\rho$ is the MPM body density, $\bm{b}$ is the body force per unit mass, $\dot{\bm{u}}$ is the velocity, $\ddot{\bm{u}}$ is the acceleration and $\bm{n}$ is the unit outward normal of the boundary $\partial\Omega_{mpm,trac}$, that represents the traction boundary. Similarly, $\partial\Omega_{mpm,disp}$ represents the portion of the boundary where displacement boundary conditions are prescribed. Constitutive equations for the Cauchy stress need to be added to complete the above set of equations. We describe the MPM solid bodies as Neo-Hookean materials, with the Cauchy stress tensor given by [28]

In Eq. (12), $\lambda$ is Lamé’s first parameter, $\bm{I}$ is the identity matrix, $\mu$ is the shear modulus, $\bm{B}=\bm{F}\,\bm{F}^{T}$ is the left Cauchy-Green strain tensor, $\bm{F}$ is the deformation gradient and $J=\det(\bm{F})$. Let $\Delta U=\{\bm{\delta u}\in\bm{H}^{1}(\Omega_{mpm})\,\,|\,\,\bm{\delta u}|_{\partial\Omega_{mpm,disp}}=0\}$. The virtual displacements are chosen as test functions for the weak formulation. The weak formulation of the momentum balance reads

where $\sigma^{s}_{ij}=\sigma_{ij}/\rho(\bm{x},t)$ is the specific Cauchy stress. The interested reader can consult [41] for more details on how to derive the above weak formulation. The integrals involved in Eq. (13) can be transformed in sums if the MPM density $\rho$ is approximated using the Dirac delta function $\delta$ as

where, $\bm{x}_{p}$ denotes the position of particle $p$, and $N_{p}$ represents the total number of particles used to discretize $\Omega_{mpm}$. Note that the number of particles $N_{p}$ is time independent, assuming that no particle is leaving $\Omega_{b}$. Substituting Eq.(14) in Eq.(13) we obtain

In the above equation, $\sigma_{ij}^{s}(\bm{x}_{p})=J(\bm{x}_{p})\,\sigma_{ij}(\bm{x}_{p})/\rho_{0}$, where $\rho_{0}$ is the density of the MPM body in the initial, undeformed configuration.

Considering the above remark, Eq. (15) can be rewritten as

The weak formulation for the MPM equations in Eq. 16 is solved numerically employing the active background grid $\mathcal{T}_{ab}$. For any given particle $p$, an interpolation from the nodes of the element $\mathcal{E}_{p}\in\mathcal{T}_{ab}$ that encloses $p$ is performed. This step is part of the assembly procedure of the background grid system and is discussed within the complete numerical algorithm.

In a time dependent framework, the position, displacement, velocity and acceleration of a given particle are functions of time. The Newmark-beta integrator is adopted for the numerical integration of $\ddot{u}_{i,p}$ in Eq. (16). If $\bm{u}^{n+1}_{p}$ denotes the particle displacement at time $t^{n+1}$ and the same notation is also adopted for particle velocity and acceleration, the method gives

for $\gamma\in[0,1]$, $\beta\in[0,0.5]$ and $\Delta t=t^{n+1}-t^{n}$. It follows that,

The complete numerical algorithm for the implicit MPM scheme presented in this paper is described next. The algorithm is divided into a series of stages and it is set in a parallel framework where different processes handle their own share of the background grid and particles.

• •

  Initialization of the MPM particles.
  The first step consists of particle initialization. We assume that the total number of particles $N_{p}$, the volume $V_{p}$ associated with each particle, and the initial density $\rho_{0}$ of the undeformed MPM body are given. With this information, the mass of each particle is computed as $m_{p}=\rho_{0}\,V_{p}$. The initial position $\bm{x}_{p}$ of each particle $p$ is also known. The values of displacement $\bm{u}_{p}$ and velocity $\dot{\bm{u}}_{p}$ of the particles are initialized as in Eq. (11). The acceleration $\ddot{\bm{u}}_{p}$ is set to zero, while the deformation gradient $\bm{F}_{p}$ at the particle is initialized to a given initial value $\bm{F}^{0}_{p}\equiv\bm{F}_{p}(0)$. Once a particle has been placed on the FEM background grid, the element $\mathcal{E}_{p}$ on the background grid that hosts the particle is determined, using the point locating method described in [11]. As already mentioned above, such an element may be different over time, so for a given $p$, $\mathcal{E}_{p}$ is actually a function of time. For ease of notation, the time dependence is not made explicit. While for each $p$ the value of $\mathcal{E}_{p}$ is known to all processes, the local coordinates $\bm{\xi}_{p}$ of $p$ are computed only by the process that owns it, so the other processes do not possess this information. To obtain the local coordinates $\bm{\xi}_{p}$ in the FEM reference configuration, an inverse mapping obtained with Newton’s method is employed. For more details on this procedure please see [11].
  

• •

  Assembly and solution of the background grid system.
  For this step, it is assumed that the particle instances at time $t^{n}$ are known and that those at time $t^{n+1}$ have to be determined. Let $\mathcal{I}_{b}=\{1,\dots,N_{b}\}$ be the set of grid points of $\mathcal{T}_{b}$. Let $\mathcal{I}_{ab}$ be the set of grid points associated with $\mathcal{T}_{ab}$ and $\mathcal{I}_{ib}=\mathcal{I}_{b}\setminus\mathcal{I}_{ab}$. The assembly procedure is carried out differently depending on whether a node of the background grid belongs to the active or inactive background grid. Both assembly procedures are described in detail in the next steps.

  1. 1.

     Assembly of the active background grid contributions.
     A loop on the particles is executed to assemble the MPM contributions. Given a particle $p$ from the loop, the element $\mathcal{E}_{p}$ on the active background grid $\mathcal{T}_{ab}$ that is currently hosting the particle is extracted. In a sense, this procedure is similar to a standard FEM procedure since the assembly is still performed in an element-by-element fashion. The differences are that now particles are used as quadrature points instead of Gauss points, and that the quadrature contributions on a given element are not computed sequentially, since particles that belong to the same element may not be all grouped together in the particle loop.

     Remark 2.

     It is important to point out that, unlike other existing implicit MPM formulations [17, 37, 32, 27], our assembly procedure does not start with a mapping from the particles to the grid nodes. The Newmark scheme in Eq.(19) is used directly on the particle instances instead of the grid instances as existing strategies do. For this reason, no mapping from the particles to the background grid is necessary.

     According to a FEM approach, the particle virtual displacement in Eq.(16) can be obtained through interpolation from the nodes of $\mathcal{E}_{p}$. If $N_{\mathcal{E}_{p}}$ denotes the number of nodes of element $\mathcal{E}_{p}$, we have $\delta u_{i,p}=\sum_{I_{\mathcal{E}_{p}}=1}^{N_{\mathcal{E}_{p}}}\delta u_{i,I_{\mathcal{E}_{p}}}\phi_{I_{\mathcal{E}_{p}}}(\bm{x}_{p})$, where $\delta u_{i,I_{\mathcal{E}_{p}}}$ represents the nodal value of the virtual displacement at node $I_{\mathcal{E}_{p}}$ and $\phi_{I_{\mathcal{E}_{p}}}$ is the finite element shape function associated with node $I_{\mathcal{E}_{p}}$. Thanks to the arbitrariness of $\delta u_{i,I_{\mathcal{E}_{p}}}$ and the fact that by definition $\delta u_{i,I_{\mathcal{E}_{p}}}=0$ on ${\mathcal{E}_{p}\cap\partial\Omega_{mpm,disp}}$, the equation for the $I_{\mathcal{E}_{p}}$-th entry of the local residual vector associated with MPM contributions is obtained from Eq. (16)

     	$\displaystyle r_{mpm,i,I_{\mathcal{E}_{p}}}$	$\displaystyle\equiv\sum_{p=1}^{N_{p}}m_{p}\Big{[}\ddot{u}_{i,p}\phi_{I_{\mathcal{E}_{p}}}(\bm{x}_{p})+\Big{(}\sum_{j=1}^{d}\sigma^{s}_{ij,p}\dfrac{\partial\phi_{I_{\mathcal{E}_{p}}}}{\partial x_{j}}(\bm{x}_{p})\Big{)}$		(20)
     		$\displaystyle-b_{i,p}\phi_{I_{\mathcal{E}_{p}}}(\bm{x}_{p})\Big{]}=0,\quad I_{\mathcal{E}_{p}}=1,\ldots N_{\mathcal{E}_{p}}.$

     In general, it may not be necessary to interpolate from the grid nodes to the particles to evaluate the body force at the particles. For instance, if the body force is a gravitational force, then no interpolation is necessary and its value at the particles can be computed directly. Therefore, the values of the body force remain evaluated at the particles, and interpolation may be performed if needed. The term $\ddot{u}_{i,p}$ in Eq. (16) is computed with the Newmark scheme in Eq. (19). Specifically, $\ddot{u}_{i,p}$ is given by

     $\displaystyle\ddot{u}_{i,p}$

     $\displaystyle=\dfrac{1}{\beta\,\Delta t^{2}}\,\widetilde{u}_{i,g}(\bm{x}_{p})-\dfrac{1}{\beta\Delta t}\,\dot{u}^{n}_{i,p}-\dfrac{1-2\beta}{2\beta}\,\ddot{u}^{n}_{i,p},$

     (21)

     where $\widetilde{u}_{i,g}(\bm{x}_{p})$ denotes the reference background grid displacement value at $\bm{x}_{p}$, defined as

     $\displaystyle\widetilde{u}_{i,g}(\bm{x}_{p})\equiv\sum_{I_{\mathcal{E}_{p}}=1}^{N_{\mathcal{E}_{p}}}\widetilde{u}_{i,I_{\mathcal{E}_{p}}}\widetilde{\phi}_{I_{\mathcal{E}_{p}}}(\bm{x}_{p}),$

     (22)

     where $\widetilde{u}_{i,I_{\mathcal{E}_{p}}}$ is as in Eq. (10). The value of $\sigma^{s}_{ij,p}$ in Eq. (20) refers to the specific Cauchy stress obtained with a partial interpolation from the grid nodes, as in [17]. Recalling the constitutive law in Eq. (12), $\sigma^{s}_{ij,p}$ is given by

     $\displaystyle\bm{\sigma}^{s}_{p}=\dfrac{\lambda\,\,\log(J_{p})}{\rho\,\,J_{p}}\,\bm{I}+\dfrac{\mu}{J_{p}}\,\Big{(}\bm{B}_{p}-\bm{I}\Big{)},$

     (23)

     where $J_{p}=\det(\bm{F}_{p})$, $\bm{B}_{p}=\bm{F}_{p}\,(\bm{F}_{p})^{T}$ and

     $\displaystyle\bm{F}_{p}=\Big{(}\widetilde{\nabla}\widetilde{\bm{u}}_{g}(\bm{x}_{p})+\bm{I}\Big{)}\,\bm{F}_{p}^{n}.$

     (24)

     In Eq. (24), $\widetilde{\nabla}$ represents the Del operator with respect to the reference configuration and $\widetilde{\nabla}\widetilde{\bm{u}}_{g}(\bm{x}_{p})$ is the background grid displacement gradient with respect to the reference configuration, evaluated at $\bm{x}_{p}$, defined as

     $\displaystyle\widetilde{\nabla}\widetilde{u}_{ij,g}(\bm{x}_{p})\equiv\sum_{I_{\mathcal{E}_{p}}=1}^{N_{\mathcal{E}_{p}}}\widetilde{u}_{i,I_{\mathcal{E}_{p}}}\dfrac{\partial\widetilde{\phi}_{I_{\mathcal{E}_{p}}}}{\partial\widetilde{x}_{j}}(\bm{x}_{p}).$

     (25)

     It is important to stress that if no coupling is considered, the undeformed configuration coincides with the reference configuration, since the reset displacement field is identically zero. This can be easily recalled considering Eq. (2).

     The next step consists of determining the scaling factors to be used in the assembly of the soft stiffness matrix. Such factors depend on flags $M_{\mathcal{E}}$ assigned to each element of the active background grid $\mathcal{T}_{ab}$. The value of a flag depends on the element’s relative position with respect to the MPM particles. Let $N_{\mathcal{E}}$ be the total number of degrees of freedom associated with an element $\mathcal{E}$ on the active background grid. For instance, for bi-quadratic quadrilateral elements $N_{\mathcal{E}}=9$. See for instance Appendix A in [11] for a sketch of the Lagrangian bi-quadratic triangle and quadrilateral with associated degrees of freedom. For simplicity, we only consider the bi-quadratic case, but the following reasoning can be easily extended to linear and serendipity elements, as well as to higher dimensions. The value of the flag $M_{\mathcal{E}}$ is initialized to $N_{\mathcal{E}}$ for all elements in the active background grid. Let $\mathcal{E}$ be any of such elements. Then the value of $M_{\mathcal{E}}$ is modified according to the following criterion: for any node $I_{\mathcal{E}}$ of $\mathcal{E}$ that also belongs to an element of $\mathcal{T}_{ib}$, $M_{\mathcal{E}}$ is decreased by a unit. Hence, for instance, bi-quadratic quadrilateral elements that are surrounded by elements that do not enclose particles will have $M_{\mathcal{E}}=1$.
     

     Next, a loop on all elements that are part of the active background grid is carried out in the standard finite element fashion, to compute the soft stiffness contribution of the local residual

     $\displaystyle{r}_{soft,i,I_{\mathcal{E}}}\equiv\mu\,C_{\mathcal{E}}\int_{\mathcal{E}}\Big{[}\widetilde{\nabla}\bm{\widetilde{u}}+\big{(}\widetilde{\nabla}\bm{\widetilde{u}}\big{)}^{T}\Big{]}_{i}\cdot\widetilde{\nabla}\widetilde{\phi}_{I_{\mathcal{E}}}dV,\,\mbox{for}\quad I_{\mathcal{E}}=1,\ldots,N_{\mathcal{E}}.$

     (26)

     The integral in Eq. (26) is approximated by means of a quadrature rule on the element $\mathcal{E}$, where Gauss points are used as quadrature points. The value of $C_{\mathcal{E}}$ is a measure of the magnitude of the soft stiffness contribution and it depends on the flags $M_{\mathcal{E}}$ previously determined. This is easily explained by the fact that the larger $M_{\mathcal{E}}$, the more elements filled with particles surround $\mathcal{E}$, and the smaller soft stiffness contribution is needed to stabilize the equations within $\mathcal{E}$. As a general rule, $C_{\mathcal{E}}$ is a non-increasing function of $M_{\mathcal{E}}$, and it is zero for any element fully surrounded by elements containing particles. The soft stiffness improves stability while decreasing accuracy. Accuracy decreases because a contribution that is not associated with the MPM problem is added to the system. As mentioned in the introduction, it is necessary to strike a balance between the need for accuracy and stability. A value of $C_{\mathcal{E}}=0$ is appropriate for elements that are fully surrounded by elements with particles but it may not be sufficient for those elements on the boundary of the MPM body. This is because they may not be filled enough with particles. It is recommended to keep $C_{\mathcal{E}}$ as small as possible, for instance in the interval $[10^{-2},10^{-8}].$ In the numerical simulations, we do not change the mass of the particles, but we test our model for increasing values of the Young’s modulus $E$. As a consequence, the stiffness of the body increases, since $C_{\mathcal{E}}$ is multiplied by $\mu$, which is proportional to $E$, while the gravitational force acting on the body remains the same. If the constant $C_{\mathcal{E}}$ remains the same, the artificial stiffness proportional to $\mu\,C_{\mathcal{E}}$ would increase. Hence, it is reasonable to decrease the value of $C_{\mathcal{E}}$ at least proportionally to the increase of the Young’s modulus. A dimensionless analysis of the equation would confirm this straightforwardly. Examples of values for $C_{\mathcal{E}}$ will be given in Section 6.

     Note that Eq. (20) and Eq. (26) only involve nodes that belong to the active background grid $\mathcal{T}_{ab}$, while the equation system at this stage is actually defined and solved for all the nodes of $\mathcal{T}_{b}$. In such a way, there is no need to construct a new background grid at each instant of time, rather part of the background grid becomes active depending on the position of $\Omega_{mpm}$, as already explained above. For all nodes that belong to $\mathcal{I}_{ib}$ (so exclusively to the inactive background grid), a fictitious equation is assembled using a lumped mass matrix, as it is shown next.
     

  2. 2.

     Assembly of the inactive background grid contributions.
     Here, a loop on all the elements of the inactive background grid is carried out in a standard finite element fashion. Let $\mathcal{E}$ be any element in $\mathcal{T}_{ib}$. Then, for any node $I_{\mathcal{E}}\in\mathcal{I}_{ib}$, the local residual is computed as the product of a lumped mass matrix with the nodal displacement vector, as follows

     $\displaystyle{r}_{mass,i,I_{\mathcal{E}}}$

     $\displaystyle\equiv\int_{\widehat{\mathcal{E}}}u_{i,I_{\mathcal{E}}}\widehat{\phi}_{I_{\mathcal{E}}}dV\approx u_{i,I_{\mathcal{E}}}\Big{(}\sum_{i_{g}=1}^{N_{\mathcal{E},g}}\widehat{\omega}_{i_{g}}\widehat{\phi}_{I_{\mathcal{E}}}(\bm{x}_{i_{g}})\Big{)},$

     (27)

     for $I_{\mathcal{E}}=1,\ldots,N_{\mathcal{E}}$, $I_{\mathcal{E}}\in\mathcal{T}_{ib}$. Recall that $\mathcal{I}_{ib}$ contains only the nodes that belong exclusively to the inactive background grid, so nodes that are shared with elements that include at least one particle are excluded from the inactive background grid contribution. In Eq. (27), for $i_{g}=1,\ldots,{N_{\mathcal{E},g}}$, $\bm{x}_{i_{g}}$ is a Gauss point of the undeformed element $\widehat{\mathcal{E}}$ associated with $\mathcal{E}$, $\widehat{\omega}_{i_{g}}$ is the quadrature weight associated with $\bm{x}_{i_{g}}$, $u_{i,I_{\mathcal{E}}}$ is $i$-th component of the nodal displacement value at the local node $I_{\mathcal{E}}$, and $\widehat{\phi}_{I_{\mathcal{E}}}$ is the shape function associated with node $I_{\mathcal{E}}$ of the undeformed element $\widehat{\mathcal{E}}$.
     

  3. 3.

     Differentiation of the residual vector.
     Let $\bm{r}_{mpm,\mathcal{E}_{p}}=[r_{mpm,1,1},r_{mpm,1,2},\ldots,r_{mpm,d,N_{\mathcal{E}_{p}}}]$ denote the local residual vector associated with the MPM contributions from element $\mathcal{E}_{p}$. The global MPM residual vector is denoted by

     $$\bm{r}_{mpm}=[r_{mpm,1,1},r_{mpm,1,2},\ldots,r_{mpm,d,N_{g}}]$$

     and is obtained in the standard finite element fashion, adding the entries of the local MPM residual vectors associated with the same global nodes. In the same way, the global soft residual vector $\bm{r}_{soft}$ and the global mass residual vector $\bm{r}_{mass}$ are constructed. It is important from an implementational point of view to ensure that the three global residual have matching dimensions and that they are all initialized to zero. Next, the global residual vector $\bm{r}$ for the background grid system is defined as

     $\displaystyle\bm{r}=\bm{r}_{mass}+\bm{r}_{soft}+\bm{r}_{mpm}.$

     (28)

     Let $\bm{u}^{k+1}=[u^{k+1}_{1,1},u^{k+1}_{1,2},\ldots,u^{k+1}_{d,N_{g}}]$ represent the vector of displacement values at the grid nodes at iteration $k+1$, then $\bm{r}(\bm{u}^{k+1})$ is a nonlinear function of $\bm{u}^{k+1}$ and to determine the values of the displacement at the grid nodes, it is necessary to solve the nonlinear equation

     $\displaystyle\bm{r}(\bm{u}^{k+1})=0.$

     (29)

     Writing $\bm{u}^{k+1}=\bm{u}^{k}+\Delta\bm{u}^{k+1}$ and expanding Eq.(29) in a Taylor series around $\bm{u}^{k}$, the following equation is obtained

     $\displaystyle 0=\bm{r}(\bm{u}^{k+1})=\bm{r}(\bm{u}^{k})+\dfrac{\partial\bm{r}(\bm{u}^{k})}{\partial\bm{u}}\Delta\bm{u}^{k+1}+O(({\Delta\bm{u}^{k+1}})^{2}).$

     (30)

     Neglecting the higher order terms, a linear system of equations is recovered

     $\displaystyle\bm{K}(\bm{u}^{k})\Delta\bm{u}^{k+1}=-\bm{r}(\bm{u}^{k}),$

     (31)

     where $\bm{K}(\bm{u}^{k})\equiv\dfrac{\partial\bm{r}(\bm{u}^{k})}{\partial\bm{u}}$ is the background grid stiffness matrix, $\Delta\bm{u}^{k+1}$ is the solution increment vector and the forcing term $-\bm{r}(\bm{u}^{k})$ is the negative residual vector associated with the solution at the previous iteration. In the numerical implementation, the matrix $\bm{K}(\bm{u}^{k})$ is obtained using automatic differentiation provided by the Adept library [21]. Hence, once the residual has been assembled, no explicit computation has to be carried out to obtain such a matrix.
     After an appropriate reordering of the grid nodes, Eq. (31) can be written as the following block system

     $\displaystyle\begin{bmatrix}\bm{K}_{ab}(\bm{u}_{ab}^{k})\quad&\bm{0}\\ \\ \bm{0}\quad&\bm{M}_{ib}\\ \end{bmatrix}\begin{bmatrix}\Delta\bm{u}_{ab}^{k+1}\\ \\ \Delta\bm{u}_{ib}^{k+1}\end{bmatrix}=\begin{bmatrix}-\bm{r}_{ab}(\bm{u}_{ab}^{k})\\ \\ -\bm{r}_{ib}(\bm{u}_{ib}^{k})\end{bmatrix},$

     (32)

     where

     		$\displaystyle\bm{K}_{ab}(\bm{u}_{ab}^{k})=\bm{K}_{soft}(\bm{u}_{ab}^{k})+\bm{K}_{mpm}(\bm{u}_{ab}^{k}),$		(33)
     		$\displaystyle\bm{r}_{ab}(\bm{u}_{ab}^{k})=\bm{r}_{soft}(\bm{u}_{ab}^{k})+\bm{r}_{mpm}(\bm{u}_{ab}^{k}),$
     		$\displaystyle\bm{K}_{soft}(\bm{u}_{ab}^{k})\equiv\dfrac{\partial\bm{r}_{soft}(\bm{u}_{ab}^{k})}{\partial\bm{u}_{ab}},\quad\bm{K}_{mpm}(\bm{u}_{ab}^{k})\equiv\dfrac{\partial\bm{r}_{mpm}(\bm{u}_{ab}^{k})}{\partial\bm{u}_{ab}},$
     		$\displaystyle\bm{r}_{ib}(\bm{u}_{ib}^{k})=\bm{r}_{mass}(\bm{u}_{ib}^{k}),\quad\bm{M}_{ib}=\dfrac{\partial\bm{r}_{mass}(\bm{u}_{ib}^{k})}{\partial\bm{u}_{ib}}.$

  4. 4.

     Solution of the background grid system.
     Let $\Delta\bm{u}^{k+1}$ denote the solution of system (32). With such an information, the displacement at the grid is advanced $\bm{u}^{k+1}=\bm{u}^{k}+\Delta\bm{u}^{k+1}$. If the displacement falls below a prescribed tolerance, the algorithm moves to the next stage, otherwise further iterations are carried out.
     

• •

  Update of Particle Instances.
  In the last stage of the algorithm, the position, velocity, acceleration and deformation gradient of the particles are updated. As for the assembly procedure, this step begins with a loop on the particles. Then, for a given particle $p$ in the loop, the element $\mathcal{E}_{p}$ that hosts $p$ is again extracted. The deformation gradient is updated in the following way

  $\displaystyle\bm{F}_{p}^{n+1}=(\widetilde{\nabla}\bm{\widetilde{u}}^{n+1}_{g}(\bm{x}_{p})+\bm{I})\,\bm{F}_{p}^{n},$

  (34)

  where

  $\displaystyle\bm{\widetilde{u}}_{g}^{n+1}(\bm{x}_{p})\equiv\sum_{I_{\mathcal{E}_{p}}=1}^{N_{\mathcal{E}_{p}}}\bm{\widetilde{u}}^{n+1}_{I_{\mathcal{E}_{p}}}\widetilde{\phi}_{I_{\mathcal{E}_{p}}}(\bm{x}_{p}).$

  (35)

  The particle displacement is obtained via interpolation from the displacement $\widetilde{\bm{u}}_{I}$ at the nodes of $\mathcal{E}_{p}$. Namely, $\bm{u}^{n+1}_{p}\equiv\bm{\widetilde{u}}_{g}^{n+1}(\bm{x}_{p})$. With this information and the particle instances at the previous instant of time, velocity and acceleration are updated using the Newmark scheme in Eq. (17) and Eq. (19) respectively. Finally, the particle positions are updated as $\bm{x}^{n+1}_{p}=\bm{x}^{n}_{p}+\bm{u}^{n+1}_{p}$.