# A stable poro-mechanical formulation for Material Point Methods leveraging
overlapping meshes and multi-field ghost penalisation

Giuliano Pretti  Robert E. Bird  Nathan D. Gavin  William M. Coombs  Charles
E. Augarde  Department of Engineering, Durham University
Science Site, South Road, Durham, DH1 3LE, UK.

###### Abstract

The Material Point Method (MPM) is widely used to analyse coupled (solid-
water) problems under large deformations/displacements. However, if not
addressed carefully, MPM u-p formulations for poro-mechanics can be affected
by two major sources of instability. Firstly, inf-sup condition violation can
arise when the spaces for the displacement and pressure fields are not chosen
correctly, resulting in an unstable pressure field. Secondly, the intrinsic
nature of particle-based discretisation makes the MPM an unfitted mesh-based
method, which can affect the system’s condition number and solvability,
particularly when background mesh elements are poorly populated. This work
proposes a solution to both problems. The inf-sup condition is avoided using
two overlapping meshes, a coarser one for the pressure and a finer one for the
displacement. This approach does not require stabilisation of the primary
equations since it is stable by design and is particularly valuable for low-
order shape functions. As for the system’s poor condition number, a face ghost
penalisation method is added to both the primary equations, which constitutes
a novelty in the context of MPM mixed formulations. This study frequently
makes use of the theories of functional analysis or the unfitted Finite
Element Method (FEM). Although these theories may not directly apply to the
MPM, they provide a robust and logical basis for the research. These
rationales are further supported by three numerical examples, which encompass
both elastic and elasto-plastic cases and drained and undrained conditions.

###### keywords:

Material Point Method , poro-mechanics , finite strain , large deformation
analyses , stability , inf-sup condition , mixed formulations

††journal: arXiv

##  1 Introduction

Since its initial publications [1, 2], the Material Point Method (MPM) has
become a widely used method for modelling solid materials undergoing extreme
deformations while maintaining a Lagrangian description of the equations. This
versatility has been tested in several engineering applications, including
snow avalanches [3, 4], ice dynamics [5, 6], slope stability [7, 8], soft
robots [9], and biomedical applications [10], to name a few. However, the MPM
has also been appreciated beyond its engineering purposes, especially in
computer graphics simulations [11, 12].

Among the different applications, fields such as geotechnics and biomechanics
deal with porous solid materials, whose mechanical behaviour is strongly
influenced by the presence of an interstitial fluid. These materials are the
subject of the study of poromechanics [13]. MPM-based analyses investigating
the behaviour of these materials have dramatically increased in the last
years, as the recent review paper by Zheng _et al._ [14] demonstrates. The
choices available to run these studies are numerous, and they can be grouped
into two categories (see, in this regards, Soga _et al._ [15]). The use of
different primary equations distinguishes between $\boldsymbol{u}-p^{(f)}$,
$\boldsymbol{v}-\boldsymbol{w}$, and $\boldsymbol{u}-p^{(f)}-\boldsymbol{w}$
formulations, where these formulations are labelled after the name of the
primary unknowns ($\boldsymbol{u}$ or $\boldsymbol{v}$ is the displacement or
velocity of the solid body, $p^{(f)}$ is the fluid pressure, and
$\boldsymbol{w}$ the fluid velocity). Since MPM formulations are based on a
point-based discretisation (these are indeed the Material Points, MPs), it is
possible to choose how many sets of MPs can be employed in a poro-mechanical
simulation. In one case, one set of MPs can be considered, simultaneously
keeping track of both phases of the mixed body. In the other, the solid and
the fluid phases are separately represented by two sets of MPs. The best
application of each combination of primary variables and the number of MP sets
firmly depends upon the considered problem, and assessing these is beyond the
scope of this work, which adopts a one-set $\boldsymbol{u}-p^{(f)}$
formulation. The reader interested in a detailed explanation and
classification of the above options can refer to these recent theses and
references therein [16, 17, 18].

One-set $\boldsymbol{u}-p^{(f)}$ formulation of the MPM can present two
sources of instabilities following the discretisation process, which can
enormously incapacitate the fundamental scope of the analyses. The first,
exhibited particularly in undrained or nearly-undrained conditions, is due to
the combination of spaces for the the displacement and the pressure fields,
and is the well-known _inf-sup_ or _LBB_ (Ladyzhenskaya-Babuška-Brezzi)
condition. After briefly introducing mixture theory in Section 2, Section 3 is
devoted to illustrate the necessary requirements to avoid this instability,
and in Section 4.4 a straightforward remedy is proposed, which is particularly
attractive for low-order interpolating functions. A MP-based discretisation
(Section 4.1) inevitably leads to small overlaps between the stencils of the
shape functions and the physical domain of the MPs (and their role as
quadrature points) and this is the source of the second instability (Section
4.2). In the literature, this problem often goes under the name of the _small-
cut_ issue. If the overlap between the shape function’s stencil and the
physical MP-based domain is practically infinitesimal, the system matrix can
become ill-conditioned and the requirements discussed in Section 3 are not
met. A possible remedy coming from the unfitted FEM literature is therefore
evaluated in Section 4.3. It must be noted that these two sources of
instabilities are not restricted to the MPM and to a $\boldsymbol{u}-p^{(f)}$
formulation for poro-mechanics: similar findings have been pointed out, for
instance, by Burman and Hansbo [19] in the context of unfitted FEM for Stokes
problems and the applicability of the methodologies discussed in this
manuscript go well beyond the MPM.

The lack of a formal analysis (see, for instance, Nguyen _et al._ [20]) of the
MPM inhibits some results, which are pivotal and well-acknowledged in the FEM
theory. However, the similarities between the two methods can be quite
substantial in some aspects. This point has been raised by different authors
(see, for instance, [21, 22]), especially in light of seeing the MPM as a FEM
where quadrature points are allowed to move independently from the mesh across
different time-steps. Exploiting these similarities in some circumstances
permits us to refer to functional analysis and the (unfitted) FEM theory and
inherit results which, even if not formally applicable, are practically
helpful in the MPM context. This idea of falling back to the FEM theory for
MPM analyses, even if not rigorous, is motivated by rationales provided
throughout the paper and supported by numerical examples in Section 5.

##  2 Fundamentals of mixture theory

To formalise the mechanics of a two-phase (solid/fluid) material modelled with
one-set of MPs, this section briefly details the underlying continuum
formulation for mixture theory.

###  2.1 Kinematics of the phases

The kinematics assumptions considered in this work are as follows:

  1. 1.

the fully saturated porous material is treated as two (solid and fluid)
juxtaposed continua;

  2. 2.

for simplicity, only one fluid material constitutes this phase;

  3. 3.

thermal effects and viscosity are neglected; and

  4. 4.

the porous material undergoes finite strains and rotations in the elasto-
plastic regime.

Having introduced the founding assumptions, let the mixed body $\mathcal{B}$
occupy a volume $\Omega$ in the original configuration (at time $t=0$) and a
volume $\omega$ in the current configuration (at the generic time $t$). Owing
to Assumption 1, each infinitesimal volume belonging to a mixed body contains
a solid and a unique (as implied by Assumption 2) fluid phase. These phases
are mapped through the different configurations via the invertible mappings
$\boldsymbol{\varphi}^{(ph)}\left(\boldsymbol{X}^{(ph)},t\right),$ with
$ph=s,f$ (standing for _s_ olid and _f_ luid), and $\boldsymbol{X}^{(ph)}$
indicating the initial position of each phase. These continua can evolve
differently in different configurations over time, which raises the question
of which phase to trace. In this sense, the mapping of the solid configuration
is _de facto_ privileged over the mapping of the fluid phase, which is mostly
undirectly traced by considering the relative velocity between the phases.
Given this consideration, the solid material particle of initial mixed volume
$d\Omega$ occupies a volume $d\omega=J\,d\Omega$ in the current configuration,
where the Jacobian $J$ is given by the determinant of the solid phase
deformation gradient111 For the sake of completeness, a second Jacobian
relative to the fluid phase could be introduced, this being given by
$J^{(f)}=\det\boldsymbol{F}^{(f)}$, with
$F^{(f)}_{iI}=\partial\varphi^{(f)}_{i}/\partial X^{(f)}_{I}$. However, since
the mixture theory assumes that the current position and mixed volume are
shared by the two phase particles, i.e.,
$\boldsymbol{x}=\varphi^{(s)}(\boldsymbol{X}^{(s)},t)=\varphi^{(f)}(\boldsymbol{X}^{(f)},t)$
and $J=J^{(f)}$, the fluid deformation gradient becomes mostly redundant. ,
this being

$\boldsymbol{F}\cdot d\boldsymbol{X}^{(s)}=d\boldsymbol{x},\quad\text{with}\
F_{iI}=\frac{\partial\varphi^{(s)}_{i}}{\partial X^{(s)}_{I}}.$ (1)

In the above equation, $\boldsymbol{x}$ denotes the material particle position
in the current configuration obtained using the solid mapping (i.e.,
$\boldsymbol{x}=\boldsymbol{\varphi}^{(s)}(\boldsymbol{X}^{(s)},t)$). The
indices $i=1,\dots,3$ and $I=1,\dots,3$ indicate the components of the
matching Cartesian basis vectors $\boldsymbol{e}_{i}=\boldsymbol{E}_{I}$ in
the current ($\boldsymbol{e}_{i}$) and reference ($\boldsymbol{E}_{I})$
configurations.

As proposed by Kröner [23], Lee [24] and Mandel [25] and considered in
Assumption 4, the deformation gradient is multiplicatively decoupled into an
elastic and plastic part as follows

$\boldsymbol{F}=\boldsymbol{F}^{e}\cdot\boldsymbol{F}^{p},$ (2)

where , due to Assumption 3, no other effects contribute to the deformation
gradient. Other measures of strain used in this work are based on the
deformation gradient (or its elastic/plastic part), such as the _left Cauchy-
Green strain_ $\boldsymbol{b}\coloneqq\boldsymbol{F}\cdot\boldsymbol{F}^{T}$
and the logarithmic strain
$\boldsymbol{\epsilon}\coloneqq\frac{1}{2}\ln\boldsymbol{b}$.

To describe the volume fractions of each phase, the Eulerian porosity $n$ is
defined as the ratio between the current fluid volume of the material particle
$d\omega^{(f)}$ and its current mixed volume $d\omega$, i.e.,
$n\coloneqq\frac{d\omega^{(f)}}{d\omega}$. _E contrario_ , $\left(1-n\right)$
gives the volume fraction of the material particle occupied by the solid
phase, i.e., $\left(1-n\right)=\frac{d\omega^{(s)}}{d\omega}$. The initial
value of the Eulerian porosity is denoted by $n_{0}$. As can be seen, the
Eulerian porosity satisfies the property

$0<n<1,$ (3)

where the extremes are excluded to always consider coexisting phases in the
material.

The _displacements_ are defined as the difference between the current and the
original positions. Given the choice of matching Cartesian basis vectors,
displacements can be defined as

$\boldsymbol{u}^{(ph)}\coloneqq\boldsymbol{\varphi}^{(ph)}\left(\boldsymbol{X}^{(ph)},t\right)-\boldsymbol{X}^{(ph)}.$
(4)

The total time derivative of the current position gives the _material
velocities_

$\boldsymbol{v}^{(ph)}\coloneqq\frac{d}{dt}\biggm{|}_{(ph)}\boldsymbol{\varphi}\left(\boldsymbol{X}^{(ph)},t\right),$
(5)

where $\frac{d}{dt}\bigm{|}_{(ph)}\left(\bullet\right)$ indicates the time
derivative of $\left(\bullet\right)$ following the $ph-$phase.

###  2.2 Mass conservation

On top of the principle of mass conservation, these further assumptions are
introduced:

  1. 5.

the constituents do not exchange mass;

  2. 6.

both constituents are incompressible; and

  3. 7.

the porosity network is homogeneous and connected across the material.

Assumption 5 permits writing the conservation of the mixture separately for
each phase, i.e.,

$\displaystyle\frac{d}{dt}\biggm{|}_{(s)}\int_{\omega^{(s)}}\rho^{(s)}\left(1-n\right)dv$
$\displaystyle=0;$ (6)
$\displaystyle\frac{d}{dt}\biggm{|}_{(f)}\int_{\omega^{(f)}}\rho^{(f)}n\,dv$
$\displaystyle=0,$ (7)

where $\rho^{(ph)}$ indicates the _mesoscopic density_ of the $ph-$phase,
which, according to Assumption 6, is constant, i.e.,
$\rho^{(ph)}=\rho_{0}^{(ph)}$. The exclusion of a double porosity network by
Assumption 7 allows the fluid mass balance to be written as in Eq. (7), i.e.,
considering only one type of porosity, approximatively homogeneous in size and
shape in the medium under investigation. Owing to the arbitrary nature of the
integration volume and the incompressibility of the solid phase (part of
Assumption 6), Eq. (6) can also be expressed as follows

$n=1-\frac{1}{J}\left(1-n_{0}\right).$ (8)

Two other equations are employed to express Eq. (7) in a convenient way. The
former is given by the relationship between material time derivatives relative
to the two different phases (see Thorpe [26]), i.e.,

$\frac{d}{dt}\biggm{|}_{(f)}\left(\bullet\right)=\frac{d}{dt}\biggm{|}_{(s)}\left(\bullet\right)+\text{{grad}}\left(\bullet\right)\cdot\left(\boldsymbol{v}^{(f)}-\boldsymbol{v}^{(s)}\right),$
(9)

where grad indicates the gradient with respect to the current position
$\boldsymbol{x}$. The second equation stems from the property relative to the
derivation of determinants of second-order tensors, and it is given by

$\frac{d}{dt}\biggm{|}_{(ph)}J=J\,\text{{div}}\cdot\boldsymbol{v}^{(ph)},$
(10)

where div is the divergence with respect to the current position
$\boldsymbol{x}$. Invoking the arbitrary nature of the integration domain
again, and using Eqs. (8)-(10), the fluid mass conservation given by Eq. (7)
can be re-written as

$\rho^{(f)}_{0}\frac{dJ}{dt}\Bigm{|}_{(s)}+J\,\text{{div}}\cdot\boldsymbol{q}^{(f)}=\rho^{(f)}_{0}\dot{J}+J\,\text{{div}}\cdot\boldsymbol{q}^{(f)}=0,$
(11)

where the shorthand notation
$\dot{\left(\bullet\right)}=\frac{d}{dt}\bigm{|}_{(s)}\left(\bullet\right)$
and $\boldsymbol{q}^{(f)}$ is the relative flux defined as

$\boldsymbol{q}^{(f)}\coloneqq\rho^{(f)}_{0}\,n\left(\boldsymbol{v}^{(f)}-\boldsymbol{v}^{(s)}\right).$
(12)

###  2.3 Balance of rate of linear momentum

Let us further assume that

  1. 8.

inertia effects are negligible.

Under this Assumption 8, the balance of the rate of the linear momentum for
the mixed body can be expressed as follows

$\text{{div}}\cdot\boldsymbol{\sigma}+\rho\ \boldsymbol{b}=\boldsymbol{0},$
(13)

where $\boldsymbol{\sigma}$ is the _total Cauchy stress_ ,
$\rho\coloneqq\rho^{(s)}_{0}\left(1-n\right)+\rho^{(f)}_{0}\,n$ is the
_current density of the porous material_ and $\boldsymbol{b}$ are the body
forces per unit mixed weight. The stress introduced in the above equation can
be further characterised. As demonstrated by Borja and Alarcón [27],
considering a perfect fluid (included in Assumption 3) and an incompressible
solid phase (Assumption 6) allows the _Terzaghi effective stress
decomposition_ 222Eq. (14) assumes the continuum mechanics sign convention
(positive pressure indicates tensile behaviour) as opposed to the geotechnical
one (positive pressure designates compressive behaviour). to hold even in the
finite strain context, this being

$\boldsymbol{\sigma}=\boldsymbol{\sigma}^{\prime}-p^{(f)}\boldsymbol{1}.$ (14)

In the above equation and throughout this manuscript, the dash
$\left(\bullet\right)^{\prime}$ denotes the _effective_ quantities of stress,
i.e., those obeying to a constitutive relationship with the the solid strains
$\boldsymbol{\sigma}^{\prime}(\boldsymbol{u}^{(s)})$, while $p^{(f)}$ is the
_Cauchy fluid pressure_.

###  2.4 Constitutive relationships

To introduce the necessary constitutive relationships, it is assumed that

  1. 9.

the material is isotropic; and

  2. 10.

the fluid flow exhibits Low Reynolds number;

Assumption 9 has two consequences. On the one hand, it results in _isotropic_
elasto-plastic behaviour relating the effective stress and the solid strains.
In particular, this work adopts the _improved Hencky material_ suggested in
Pretti _et al._ [28] to avoid negative values of the Eulerian porosity. The
effective free energy function $\Psi^{\prime}$ of this material is given by

$\Psi^{\prime}\left(\boldsymbol{\epsilon},\boldsymbol{\alpha}\right)=\frac{K}{2\,n}\left(\epsilon_{v}^{e}\right)^{2}+\frac{3}{2}G\left(\epsilon_{q}^{e}\right)^{2}+\tilde{\Psi}^{\prime}\left(\boldsymbol{\alpha}\right),$
(15)

where $K>0$ and $G>0$ are the bulk parameter and the shear modulus, while
$\tilde{\Psi}^{\prime}$ defines a part of the free energy function responsible
for the kinematic hardening, based on a set of internal variables
$\boldsymbol{\alpha}$. A few invariants of the logarithmic strains have also
been employed in the above equation, which are defined as

$\epsilon_{v}\coloneqq\boldsymbol{\epsilon}\boldsymbol{:}\boldsymbol{1};\qquad\boldsymbol{e}\coloneqq\boldsymbol{\epsilon}-\frac{\epsilon_{v}}{3}\boldsymbol{1};\qquad\epsilon_{q}\coloneqq\sqrt{\frac{2}{3}\,\boldsymbol{e}\boldsymbol{:}\boldsymbol{e}},$
(16)

with $\boldsymbol{:}$ being the double contraction operator between tensors.
If isotropic permeability is assumed (second consequence of Assumption 7)
together with Assumption 10, the _Darcy-Weisbach_ law relates the relative
flux to the fluid pressure gradient and body forces as follows

$\boldsymbol{q}^{(f)}=-\frac{\upkappa}{g}\left(\text{{grad}}\,p^{(f)}-\rho^{(f)}\boldsymbol{b}\right),$
(17)

with $\upkappa$ being the _hydraulic conductivity_ and $g$ the magnitude of
the gravitational force. The _Kozeny-Carman_ formula is assumed to account for
the dependency of the hydraulic conductivity from the porosity, i.e.,

$\upkappa=c_{1}\frac{n^{3}}{\left(1-n\right)^{2}},$ (18)

with $c_{1}$ being a constant parameter. As embedded in Assumption 6, the
fluid pressure is unrelated to any constitutive relationship and, as pointed
out by Miehe _et al._ [29], acts as a Lagrange multiplier by shifting the
dependency of stresses on strain (i.e., the idea of effective stress)
depending on the water mass conservation constraint (i.e., drained/undrained
conditions). This observation has direct consequences on the mapping processes
detailed in Section 4.6, since the fluid pressure does not relative to any
form of energy.

##  3 Continuum weak form and arising requirements

This section introduces a weak statement of the equations for a
$\boldsymbol{u}-p^{(f)}$ formulation, which results in a saddle point
formulation (Section 3.1). As is well-acknowledged, this kind of problem can
be solved only under specific conditions.

###  3.1 From strong to weak form of equations

The primary equations employed in a $\boldsymbol{u}-p^{(f)}$ formulation are
given by the mixture equilibrium Eq. (13) and the fluid mass conservation Eq.
(11), these being

$\displaystyle\text{{div}}\cdot\left(\boldsymbol{\sigma}^{\prime}-p^{(f)}\,\boldsymbol{1}\right)+\rho\
\boldsymbol{b}$ $\displaystyle=\boldsymbol{0}\quad\text{in}\ \ \omega;$ (19)
$\displaystyle\rho_{0}^{(f)}\dot{J}-J\,\text{{div}}\cdot\left(\frac{\upkappa}{g}\left(\text{{grad}}\,p^{(f)}-\rho^{(f)}\boldsymbol{b}\right)\right)$
$\displaystyle=0\quad\text{in}\ \ \omega^{(f)},$ (20)

where the Terzaghi effective stress decomposition Eq. (14) and the Darcy-
Weisbach law Eq. (17) have been substituted. The current boundary
$\gamma=\partial\omega$ of the considered mixed body is partitioned as follows

$\displaystyle\gamma=\gamma^{\bar{u}}\cup\gamma^{\bar{t}}=\gamma^{\bar{p}}\cup\gamma^{\bar{q}};$
(21)
$\displaystyle\gamma^{\bar{u}}\cap\gamma^{\bar{t}}=\emptyset=\gamma^{\bar{p}}\cap\gamma^{\bar{q}},$
(22)

where $\gamma^{(\bullet)}$ is a particular portion of the boundary. Prescribed
boundary conditions (BCs) are given on these parts and are as follows

$\displaystyle\boldsymbol{u}$ $\displaystyle=\bar{\boldsymbol{u}}\qquad$
$\displaystyle\hbox{on}\ \gamma^{\bar{u}};$ (23)
$\displaystyle\boldsymbol{\sigma}\,\boldsymbol{n}$
$\displaystyle=\bar{\boldsymbol{t}}\qquad$ $\displaystyle\hbox{on}\
\gamma^{\bar{t}};$ (24) $\displaystyle p^{(f)}$ $\displaystyle=\bar{p}\qquad$
$\displaystyle\hbox{on}\ \gamma^{\bar{p}};$ (25)
$\displaystyle\boldsymbol{q}\cdot\,\boldsymbol{n}$
$\displaystyle=\bar{q}\qquad$ $\displaystyle\hbox{on}\ \gamma^{\bar{q}}.$ (26)

The above equations constitute the Dirichlet, Eqs. (23) and (25), and the
Neumann, Eqs. (24) and (26), BCs for the mixture equilibrium (19) and fluid
mass conservation (20), respectively.

The weak form of the above problem is obtained by introducing test functions
belonging to the required function spaces and integrating over the relative
domains333For the fluid mass conservation, it must be noted that, technically,
the correct integration domain is the fluid. However, this domain and that of
the mixture are related by the Eulerian porosity, i.e.,
$d\omega^{(f)}=n\,d\omega$. If this (i.e., the porosity) is simplified for all
the terms in the fluid mass conservation, the integration over the mixture
volume is _de facto_ achieved.. By applying the divergence theorem and using
Neumann BCs (24) and (26), the weak form can be stated as follows: find
$\boldsymbol{u}\in\mathscr{V}_{\bar{\boldsymbol{u}}}$ and
$p^{(f)}\in{\mathscr{W}_{\bar{p}}}$ such that, for $t\in[0,\tilde{t}]$,

$\delta\Pi^{eq}\left(\left(\boldsymbol{u};p^{(f)}\right);\delta\boldsymbol{w}\right)\coloneqq\int_{\omega}\text{{grad}}\left(\delta\boldsymbol{w}\right)\boldsymbol{:}\left(\boldsymbol{\sigma}^{\prime}-p^{(f)}\,\boldsymbol{1}\right)\
dv-\int_{\omega}\rho\,\delta\boldsymbol{w}\cdot\boldsymbol{b}\ dv\\\
-\int_{\gamma^{\bar{t}}}\delta\boldsymbol{w}\cdot\bar{\boldsymbol{t}}\
da=0,\qquad\forall\delta\boldsymbol{w}\in\mathscr{V}_{\boldsymbol{0}};$ (27)
$\delta\Pi^{cons}\left(\left(\boldsymbol{u};p^{(f)}\right);\delta\eta\right)\coloneqq\int_{\omega}\rho_{0}^{(f)}\,\delta\eta\frac{\dot{J}}{J}\
dv+\int_{\omega}\frac{\upkappa}{g}\,\text{{grad}}\left(\delta\eta\right)\cdot\left(\text{{grad}}\,p^{(f)}-\rho^{(f)}\boldsymbol{b}\right)dv\\\
-\int_{\gamma^{\bar{q}}}\delta\eta\,\bar{q}\
da=0,\qquad\forall\delta\eta\in\mathscr{W}_{0},$ (28)

giving the initial condition
$\boldsymbol{u}\left(t=0\right)=\boldsymbol{u}_{0}$. The potentials $\Pi^{eq}$
and $\Pi^{cons}$ and their first variation
$\delta\left(\bullet\right)\left(\left(\dots\right);\delta\left(\tilde{\bullet}\right)\right)$
with respect to their argument $\left(\tilde{\bullet}\right)$ have been above
employed, while the spaces of trial functions are defined as follows

$\displaystyle\mathscr{V}_{\bar{\boldsymbol{u}}}$
$\displaystyle=\left\\{\boldsymbol{u}\in\left[H^{1}\left(\omega\right)\right]^{n^{dim}}\bigm{|}\boldsymbol{u}=\bar{\boldsymbol{u}}\
\hbox{on}\ \gamma^{\bar{u}}\right\\};$ (29)
$\displaystyle\mathscr{W}_{\bar{p}}$
$\displaystyle=\left\\{p^{(f)}\in\left[H^{1}\left(\omega\right)\right]\bigm{|}p^{(f)}=\bar{p}\
\hbox{on}\ \gamma^{\bar{p}}\right\\},$ (30)

where $H^{1}\left(\omega\right)$ denotes the Sobolev space of degree one on
$\omega$. The spaces of the test functions $\delta\boldsymbol{w}$ and
$\delta\eta$ are denoted by $\mathscr{V}_{\boldsymbol{0}}$ and
$\mathscr{W}_{0}$ and follow similar definitions as the above, with the
exception of being zero on the boundary. It is well-known (see, for instance,
Dortdivanlioglu _et al._ [30]) that the above weak problem Eqs. (27)-(28) can
be given by the stationarity of the following saddle point problem

$\inf_{\boldsymbol{w}\in\mathscr{V}_{\boldsymbol{0}}}\sup_{q\in\mathscr{W}_{0}}\Pi\left(\boldsymbol{w};\eta\right),$
(31)

with
$\Pi\left(\boldsymbol{w};\eta\right)=\Pi^{eq}\left(\boldsymbol{w};\eta\right)-\Pi^{cons}\left(\boldsymbol{w};\eta\right)$.

###  3.2 Linearisation

Eqs. (27) and (28) represent a non-linear system of equations in
$\boldsymbol{u}$. This is primarily due to having set the problem in the
context of finite strain mechanics. On top of this, other non-linearities,
such as elasto-plastic behaviour or constitutive equation (15), can be
included too. To solve these equations, the Newton-Raphson (NR) method
stipulates that linearisation of Eqs. (27) and (28) at a solution
$\bigl{(}\check{\boldsymbol{u}};\check{p}^{(f)}\bigr{)}$ is required, this
being

$\displaystyle
0=\delta\Pi\Bigl{(}\bigl{(}\check{\boldsymbol{u}};\check{p}^{(f)}\bigr{)};\left(\delta\boldsymbol{w},\delta\boldsymbol{\eta}\right)\Bigr{)}$
$\displaystyle\approx\delta\Pi\left(\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\left(\delta\boldsymbol{w},\delta\boldsymbol{\eta}\right)\right)$
$\displaystyle\qquad+\delta\Biggl{(}\biggl{(}\delta\Pi\Bigl{(}\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\left(\delta\boldsymbol{w},\delta\boldsymbol{\eta}\right)\Bigr{)}\biggr{)};\bigl{(}\delta\boldsymbol{u},\delta\boldsymbol{p}^{(f)}\bigr{)}\Biggr{)}$
$\displaystyle=\delta\Pi^{eq}\Bigl{(}\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\left(\delta\boldsymbol{w}\right)\Bigr{)}-\delta\Pi^{cons}\Bigl{(}\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\bigl{(}\delta\eta\bigr{)}\Bigr{)}$
$\displaystyle\qquad+\delta\Biggl{(}\biggl{(}\delta\Pi^{eq}\Bigl{(}\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\delta\boldsymbol{w}\Bigr{)}\biggr{)};\bigl{(}\delta\boldsymbol{u},\delta\boldsymbol{p}^{(f)}\bigr{)}\Biggr{)}$
$\displaystyle\qquad\qquad-\delta\Biggl{(}\biggl{(}\delta\Pi^{cons}\Bigl{(}\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\delta\eta\Bigr{)}\biggr{)};\bigl{(}\delta\boldsymbol{u},\delta\boldsymbol{p}^{(f)}\bigr{)}\Biggr{)}.$
(32)

The second variations appearing in the above equations can expressed as
follows

$\delta\Biggl{(}\biggl{(}\delta\Pi^{eq}\Bigl{(}\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\delta\boldsymbol{w}\Bigr{)}\biggr{)};\bigl{(}\delta\boldsymbol{u},\delta\boldsymbol{p}^{(f)}\bigr{)}\Biggr{)}-\delta\Biggl{(}\biggl{(}\delta\Pi^{cons}\Bigl{(}\bigl{(}\boldsymbol{u};p^{(f)}\bigr{)};\delta\eta\Bigr{)}\biggr{)};\bigl{(}\delta\boldsymbol{u},\delta\boldsymbol{p}^{(f)}\bigr{)}\Biggr{)}\\\
=\begin{bmatrix}\delta\boldsymbol{w},\delta\eta\end{bmatrix}\underbrace{\begin{bmatrix}\mathbf{A}&\mathbf{B}^{(1)}\\\
\mathbf{B}^{(2)}&\mathbf{C}\end{bmatrix}}_{\coloneqq\mathbf{J}}\begin{bmatrix}\delta\boldsymbol{u}\\\
\delta p^{(f)}\end{bmatrix},$ (33)

where $\mathbf{J}$ indicates the Jacobian matrix. Under the assumptions made
so far, the submatrix $\mathbf{C}$ is symmetric, i.e.,
$\mathbf{C}=\mathbf{C}^{T}.$

In the small strain regime, the above linear system satisfies also the useful
properties $\mathbf{A}=\mathbf{A}^{T}$ and
$\left(\mathbf{B}^{(1)}\right)^{T}=\mathbf{B}^{(2)}$. In this context, the
Jacobian matrix can exhibit two situations which are well-acknowledged for
saddle-point problems. On the one hand, in the case of drained or partially
drained processes, the above Jacobian is as follows

$\mathbf{J}=\begin{bmatrix}\mathbf{A}&\mathbf{B}^{T}\\\
\mathbf{B}&\mathbf{C}\end{bmatrix}.$ (34)

The solvability (uniqueness of solution) of the above linear system is
guaranteed if the bilinear forms associated with the submatrices $\mathbf{A}$
and $\mathbf{B}$ are continuous, and $\mathbf{A}$ and $\mathbf{C}$ are
coercive (see, for instance, Boffi _et al._ [31], Proposition 4.3.1). On the
other hand, the process can take place in nearly-undrained conditions.
Physically speaking, these circumstances can occur when load rates are rapidly
applied with low values of hydraulic conductivity or when BCs do not allow for
drainage. Theses cases result in the following Jacobian

$\mathbf{J}=\begin{bmatrix}\mathbf{A}&\mathbf{B}^{T}\\\
\mathbf{B}&\epsilon\,\tilde{\mathbf{C}}\end{bmatrix},$ (35)

where $\mathbf{C}=\epsilon\,\tilde{\mathbf{C}}$ with $0<\epsilon<<1$ indicates
that the submatrix $\mathbf{C}$ contains small entries compared to the other
submatrices. In the above condition, solvability is guaranteed (see, e.g.,
Boffi _et al._ [31], Theorem 4.3.4) if the bilinear form associated with
$\mathbf{A}$ and $\mathbf{C}$ are continuous and $\mathbf{A}$ is coercive.
Furthermore, the solution to Eq. (35) is stable if $\mathbf{B}^{T}$ is such
that its bilinear form satisfies the _inf-sup_ condition, i.e.,

$\inf_{\left({\delta\eta}\in{\mathscr{W}_{0}}\right)}\
\sup_{\left({\delta\boldsymbol{w}}\in{\mathscr{V}_{\boldsymbol{0}}}\right)}\dfrac{\int_{\omega}{\delta\eta}\
\text{{div}}\cdot{\delta\boldsymbol{w}}\,dv}{||{\delta\eta}||_{H^{0}}\
||{\delta\boldsymbol{w}}||_{{\mathscr{V}_{\boldsymbol{0}}}}}\geq\alpha>0,$
(36)

where the L-2 norm $||\left(\bullet\right)||_{H^{0}}$ must be considered even
for poro-mechanical problems (see the discussions in Mira _et al._ [32] and
Dortdivanlioglu _et al._ [30]). Furthermore, the use of this norm in Eq. (36)
allows use of numerical test proposed by Chapelle and Bathe [33] for
incompressible elasticity or Stokes problems.

From the above discussions on drained and undrained cases, it emerges that
three requirements are necessary for inverting the Jacobian matrix without
incurring in unstable solutions (neglecting for the moment the continuity of
bilinear forms of submatrices $\mathbf{A}$ and $\mathbf{C}$). These conditions
are the coercivity of the bilinear forms associated with $\mathbf{A}$ and
$\mathbf{C}$, and the satisfaction of the inf-sup condition for the bilinear
form associated with $\mathbf{B}$.

It must also be stressed that adopting the improved Hencky material defined by
Eq. (15) has important consequences for the coercivity of $\mathbf{C}$, since
it avoids negative values of the Eulerian porosity, as detailed in Pretti _et
al._ [28]. While this criterion does not come from the numerical requirements
of the problem under investigation, its physically based nature substantially
affects numerical outcomes. This constraint and its relative consequences hold
even in the case of compressible material. In this case, the reader is invited
to refer to the work of Nedjar [34, 35] on circumventing negative values in
Eulerian porosity. Owing to the adoption of the improved Hencky material in
this work, this physical constraint will not be discussed further.

##  4 Stable MPM discretisation

This section introduces an MPM $\boldsymbol{u}-p^{(f)}$ formulation and
discusses when the conditions for solvability discussed in the previous
section are not fulfilled. Remedies to mitigate the effects of these losses
are also proposed.

###  4.1 MPM algorithm

Figure 1: MPM step phases.

Figure 1 shows the phases in the overall MPM algorithm. The MPM-based
discretisation stipulates that the considered body of global volume $\omega$
is discretised by a cloud of Material Points (MPs), i.e.,
$\omega\approx\cup_{mp=1}^{N^{mp}}\ {{}^{mp}\omega}\coloneqq{{}^{MP}\omega}$;
the MPs carry all the information necessary to run the analysis. These points
are immersed into a (usually Cartesian) grid/mesh $\mathcal{T}$ discretising a
portion of the Euclidean space $\mathcal{E}$ of $n^{dim}$ dimensions, which
fully contains the considered body $\omega$ (Phase $(\mathrm{i})$ in Figure
1).

Finite-dimensional test and trial functions are defined on a mesh of generic
length $\tilde{h}$. $\tilde{h}$, when employed as a superscript, denotes
finite-dimensional quantities. The trial and test functions functions are
interpolated with the help of nodal shape functions, similar to the FEM. The
finite-dimensional spaces of the trial functions are given by

${{}^{\tilde{h}}\mathscr{V}_{\bar{\boldsymbol{u}}}}=\Bigl{\\{}{{}^{\tilde{h}}\boldsymbol{u}}\in\left[C^{m}\left({}^{\tilde{h}}\bar{\omega}\right)\right]^{n^{dim}}\bigm{|}\,{{}^{\tilde{h}}\boldsymbol{u}}\left(\boldsymbol{x}\right)=\mathrm{N}^{\tilde{h},u}_{a}\left(\boldsymbol{x}\right)\mathbf{u}_{a}\
\forall a\in\text{clos}\left({{}^{\tilde{h}}\mathcal{T}^{act}}\right),\,\\\
{{}^{\tilde{h}}\boldsymbol{u}}=\bar{\boldsymbol{u}}\ \hbox{on}\
{{}^{\tilde{h}}\bar{\gamma}^{\bar{u}}}\Bigr{\\}};$ (37)
${{}^{\tilde{h}}\mathscr{W}_{\bar{p}}}=\left\\{{}^{\tilde{h}}p^{(f)}\in
C^{m}\left({}^{\tilde{h}}\bar{\omega}\right)\bigm{|}\,{{}^{\tilde{h}}p^{(f)}}\left(\boldsymbol{x}\right)=\mathrm{N}^{\tilde{h},p}_{a}\left(\boldsymbol{x}\right)\,\mathrm{p}_{a}\
\forall
a\in\text{clos}\left({{}^{\tilde{h}}\mathcal{T}^{act}}\right)\right\\},$ (38)

where
${}^{\tilde{h}}\bar{\omega}=\left(\bigcup_{T\in{{}^{\tilde{h}}\mathcal{T}^{act}}}T\right)$
denotes the volume defined by the set of active grid elements
${{}^{\tilde{h}}\mathcal{T}^{act}}$, i.e.,

${{}^{\tilde{h}}\mathcal{T}^{act}}=\bigl{\\{}T\in{{}^{\tilde{h}}\mathcal{T}}\bigm{|}\exists\,mp:\mathrm{N}^{\tilde{h}}_{a}\left({}^{mp}\boldsymbol{x}\right)>0\
\text{with}\,a\in\text{clos}(T)\bigr{\\}},$ (39)

with $\mathrm{N}^{\tilde{h}}_{a}$ being the generic low-order shape functions
at the $a-$th node belonging to the element $T$ and
$\text{clos}\left(\bullet\right)$ denotes the closure of the entity
$\left(\bullet\right)$444 This work often refers to the closure of a generic
entity, which is a shorthand for indicating the set of entities contained by
the considered generic one and the entity itself. If, for instance, a bi-
dimensional mesh element $T$ is considered, the closure of $T$ includes the
faces and nodes belonging to $T$ and $T$ itself. . As detailed below in
Section 4.4, this work assumes two kinds of shape functions. The first choice
is given by first-order Lagrange polynomials, and results in
$C^{0}({{}^{\tilde{h}}\bar{\omega}})$ functions (i.e., $m=0$). Since they were
first adopted in the original MPM formulation (see Sulsky _et al._ [1, 2]),
these shape functions are labelled sMPM (standard MPM). On the other hand,
GIMPM shape functions (Generalised Interpolation MPM, see, for instance, [36,
37]), resulting in piece-wise first-order and second-order polynomials, are
also considered, giving $C^{1}({{}^{\tilde{h}}\bar{\omega}})$ functions (i.e.,
$m=1$).

To initialise the grid with the required information, a mapping from the MPs
to the grid nodes is usually required (Phase $(\mathrm{ii})$ in Figure 1).
This procedure is named Point-2-Grid (P2G) mapping and, for the current
formulation, is addressed in Section 4.6.

The (standard) assembly process of the Lagrangian equations (Phase
$(\mathrm{iii})$ in Figure 1) takes place at the grid nodes and employs the
MPs as quadrature points. In the case of the considered
$\boldsymbol{u}-p^{(f)}$ formulation, the equations become (dropping the
dependencies from the position)

$\int_{{{}^{MP}\omega}}\updelta\mathbf{w}_{a}\,\text{{grad}}\left(\mathrm{N}^{\tilde{h},u}_{a}\right)\boldsymbol{:}\left(\boldsymbol{\sigma}^{\prime}-\mathrm{N}^{\tilde{h},p}_{b}\,\mathrm{p}_{b}\,\boldsymbol{1}\right)\
dv-\int_{{{}^{MP}\omega}}\left(\updelta\mathbf{w}_{a}\,\mathrm{N}^{\tilde{h},u}_{a}\right)\cdot\rho\boldsymbol{b}\
dv\\\
-\int_{\gamma^{\bar{t}}}\left(\updelta\mathbf{w}_{a}\,\mathrm{N}^{\tilde{h},u}_{a}\left(\boldsymbol{x}\right)\right)\cdot\bar{\boldsymbol{t}}\
da=0,\quad\forall\updelta\mathbf{w}_{a}\in\
^{\tilde{h}}\mathscr{V}_{\boldsymbol{0}};$ (40)
$\int_{{{}^{MP}\omega}}\left(\updelta\upeta_{a}\,\mathrm{N}^{\tilde{h},p}_{a}\right)\rho_{0}^{(f)}\frac{\dot{J}}{J}\
dv+\int_{{{}^{MP}\omega}}\frac{\upkappa}{g}\,\updelta\upeta_{a}\,\text{{grad}}\left(\mathrm{N}^{\tilde{h},p}_{a}\right)\cdot\left(\text{{grad}}\left(\mathrm{N}^{\tilde{h},p}_{b}\right)\mathrm{p}_{b}-\rho^{(f)}\boldsymbol{b}\right)dv\\\
-\int_{\gamma^{\bar{q}}}\left(\updelta\upeta_{a}\,\mathrm{N}^{\tilde{h},p}_{a}\right)\,\bar{q}\
da+\gamma^{pen}\int_{\gamma^{\bar{p}}}\left(\updelta\upeta_{a}\,\mathrm{N}^{\tilde{h},p}_{a}\right)\,\left(\mathrm{N}^{\tilde{h},p}_{b}\,\mathrm{p}_{b}-\bar{p}\right)da=0,\quad\forall\updelta\upeta_{a}\in\
{{}^{\tilde{h}}\mathscr{W}_{0}},$ (41)

where the volume integrals are approximated as follows

$\int_{{{}^{MP}\omega}}\left(\bullet\right)dv\approx\sum_{mp=1}^{N^{mp}}{{}^{mp}\omega}\left(\bullet\left({}^{mp}\boldsymbol{x}\right)\right).$
(42)

As for the boundary integrals, this work follows the procedure explained by
Bird _et al._ [38], leveraging the MPs’ corners (see Bird _et al._ [38] for
more details). This method is also exploited for non-conforming pressure
Dirichlet BCs, represented by the penalty term in Eq. (41), with
$\gamma^{pen}$ being a penalty parameter. Conversely, this work always
considers mesh-conforming displacement BCs, as assumed in Eq. (37). Eqs. (40)
and (41) are solved for the nodal unknowns (Phase $(\mathrm{iv})$ in Figure
1), and this work considers an interactive Newton-Raphson (NR) scheme based on
the discrete counterpart of Eq. (3.2). Since the updated solution lies on the
grid, a subsequent mapping from this to the MPs is required to update the
solution at the MPs level (Phase $(\mathrm{v})$ in Figure 1). This is called
the Grid-2-Point (G2P) mapping and is discussed in Section 4.6. The grid is
discarded at the end of the step (Phase $(\mathrm{vi})$ in Figure 1), and a
further step can start with an undistorted grid.

###  4.2 The small-cut issue in the MPM

As will be clear from the description of the phases of the MPM, the immersion
of the clouds of MPs into a mesh and their use as quadrature points generates
an intrinsically unfitted method since
${{}^{MP}\omega}\subseteq{{}^{\tilde{h}}\bar{\omega}}$. Moreover, an extremely
small overlap between the shape function’s stencil and the MP’s domains can
lead to ill-conditioned systems. This problem, named the small-cut issue is
exacerbated especially for those elements on the domain’s boundaries. As
demonstrated by Burman [39] in the context of unfitted FEM, this issue
provokes the loss of the coercivity for the bilinear forms computed on these
elements and the consequent impossibility of solving the relative linear (or
linearised) system of equations.

To manage this issue, Burman [39] has proposed adding a penalty term, named
the ghost penalty stabilisation, which extends the coercivity of the bilinear
forms to the partially filled boundary elements . While strictly enforcing
coercivity in an MPM formulation is difficult to prove (properties such as
continuity or coercivity rely on norms which, in turn, are based on quadrature
rules more precise than Eq. (42)), it is nonetheless possible to limit the
condition number of the matrices by the use of the ghost penalty method (see
Coombs [40]). The addition of the ghost penalty restores to some degree
coercivity on those elements suffering from the small-cut issue. Hence, while
formalisms must be dropped for the above-mentioned causes, MPM-based discrete
linear systems can be inverted independently from the small-cut issue, as long
as the ghost penalty (or a similar technique, see the discussion in Section
4.3) is introduced.

Among the ghost penalties introduced in the literature for the unfitted FEM
(see Burman _et al._ [41]), the _face_ ghost penalty555 While the _face_ ghost
penalty is the term usually adopted in the literature, this work refers to
_facet_ to indicate, in the $n^{dim}-$th dimensional space, the face with
dimension $n^{dim}-1$, i.e., a point, an edge, and a face for $n^{dim}=1,2,3$,
respectively.  has been introduced to the MPM by Coombs [40], and an extension
to a $\boldsymbol{u}-p^{(f)}$ formulation for the unfitted FEM has been
proposed by Liu _et al._ [42]. The following section extends the face ghost
stabilisation for the $\boldsymbol{u}-p^{(f)}$ MPM formulation introduced
above.

###  4.3 Face ghost stabilisation

In the case of low-order shape functions, the face ghost stabilisation states
that the following bilinear form must be added to the considered primary
equations, i.e.,

$j^{F}\left({{}^{\tilde{h}}\delta\boldsymbol{v}},{{}^{\tilde{h}}\boldsymbol{v}}\right)\coloneqq\tilde{h}\int_{{}^{\tilde{h}}\gamma^{F}}\left[\left[\left(\text{{grad}}\left({{}^{\tilde{h}}\delta\boldsymbol{v}}\right)\cdot\mathbf{n}\right)\right]\right]\cdot\left[\left[\left(\text{{grad}}\left({}^{\tilde{h}}\boldsymbol{v}\right)\cdot\mathbf{n}\right)\right]\right]\
da,$ (43)

where ${}^{\tilde{h}}\boldsymbol{v}$ and ${}^{\tilde{h}}\delta\boldsymbol{v}$
are the generic trial and test functions and
$[[\text{{grad}}\left(\bullet\right)\,\mathbf{n}]]$ is the jump in the normal
gradient across a shared facet $F$. The procedure to select these facets
(listed in $\mathcal{F}$) describing the surface ${}^{\tilde{h}}\gamma^{F}$ is
outlined in Algorithms 1 and 2. In particular, since the MPM retains no
explicit representation of the boundary, Algorithm 1 is necessary to track
which grid elements constitute the domain’s boundary
${{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}$. To read
this algorithm, it is necessary also to introduce the set of inactive elements
${}^{\tilde{h}}\mathcal{T}^{int}$, defined as
${}^{\tilde{h}}\mathcal{T}^{in}\coloneqq{{}^{\tilde{h}}\mathcal{T}}\setminus{{}^{\tilde{h}}\mathcal{T}^{act}}$.
Among the facets belonging to the set of boundary elements
${{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}$,
Algorithm 2 selects those satisfying one of the following criteria:

  * •

the facet is shared between one boundary element and one active non-boundary
element; or

  * •

the facet is between two adjacent boundary elements.

While selecting the former set of facets (criterion at the first bullet point)
is justified by prolonging the coercivity from the bulk of the material, the
latter (criterion at the second bullet point) is mainly required for boundary
corner elements not directly connected to an active non-boundary element. The
union of the selected facets, i.e.,

${}^{\tilde{h}}\gamma^{F}=\bigcup_{F\in\mathcal{F}}F$ (44)

constitutes the surface ${}^{\tilde{h}}\gamma^{F}$ necessary for Eq. (43).
Unlike the volume integrals in Eq. (42), Eq. (43) is integrated via Gauss-
Legendre quadrature rule using Gauss Points seeded on the grid.

${{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}\leftarrow[]$;

el_count $\leftarrow$ 1;

foreach  _element $T\in{{}^{\tilde{h}}\mathcal{T}^{act}}$ _ do // loop over
active elements

face_count $\leftarrow$ 1;

foreach  _facet $F\in\text{clos}(T)$ _ do // loop over facets of the selected
active element

if  _ $F\in\text{clos}({{}^{\tilde{h}}\mathcal{T}^{in}})$ _ then // check if
facet is also shared with an inactive element

Add $T$ to
${{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}$;

end if

face_count ++

end foreach

el_count ++

end foreach

Algorithm 1 Selection of the boundary elements belonging to
${{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}$.

$\mathcal{F}\leftarrow[]$;

el_count $\leftarrow$ 1;

foreach  _element
$T_{i}\in{{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}$ _
do // loop over boundary elements

face_count $\leftarrow$ 1;

foreach  _facet $F\in\text{clos}(T_{i})$ _ do // loop over facets of the
selected boundary element

if  _
$F\in\text{clos}\left({{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}\setminus
T_{i}\right)$ _ then // check if facet is shared with another boundary element

Add $F$ to $\mathcal{F}$;

end if

if  _
$F\in\text{clos}\left({{}^{\tilde{h}}\mathcal{T}^{act}}\setminus{{}^{\tilde{h}}\mathcal{T}^{\partial{{}^{\tilde{h}}}\bar{\omega}}}\right)$
_ then // check if facet is shared with active non-boundary element

Add $F$ to $\mathcal{F}$;

end if

face_count ++

end foreach

el_count ++

end foreach

Algorithm 2 Selection of the facets for the face ghost penalty.

Algorithms 1 and 2 (as well as Eq. (43)) are relevant to sMPM and GIMPM shape
functions. Unlike the unfitted FEM, in which the boundary elements constitute
a single layer of elements, high-order MPMs can contain multiple layers of
boundary elements due to the stencil of these high-order shape functions,
which cover multiple elements simultaneously. This situation can also occur
for the GIMPM. However, since only two layers of boundary elements can be
active simultaneously in the worst-case scenario for the GIMPM, Algorithms 1
and 2 prove to be sufficient (see Example 5.2).

As stabilising multiple layers can be quite burdensome, Yamaguchi _et al._
[43] adapted the Extended B-splines (EBS, see Höllig _et al._ [44, 45])
technique to stabilise MPM formulations employing B-spline shape functions.
Another interesting technique has been proposed by Badia _et al._ [46] for the
unfitted FEM, which aggregates elements on the boundary suffering from the
small-cut issue. However, these techniques can be problematic to integrate
into the current framework: on the one hand, the EBS-MPM cannot be used for
the considered low-order shape functions, while aggregating can be
particularly difficult, especially in light of the inf-sup stable overlapping
meshes adopted in this manuscript. These are explained in the following
section.

###  4.4 QkSD-Qk spatial discretisation

To comply with the inf-sup requirement, the concept of QkSD-Qk elements (where
the subscript $SD$ stands for _S_ ub-_D_ ivided, while $k$ is the degree of
the shape functions polynomials) from the FEM is here adapted to the MPM. This
is achieved by overlapping a coarser grid for the pressure field (of length
$h$) on a finer mesh for the displacement field (of length $h/2$). For
simplicity, both these meshes are considered as Cartesian.

QkSD-Qk elements are well-acknowledged in the FEM (see, for instance,
Dortdivanlioglu _et al._ [30] for a $\boldsymbol{u}-p^{f}$ formulation), and
the idea can be traced at least back to Bercovier and Pironneau [47], who gave
the first error-estimate for two- and three-dimensional FEM of these kind of
elements. In the case of a MPM formulation, to the best of authors’ knowledge,
QkSD-Qk concept has been applied only to problems involving incompressible
elasticity and finite elements with B-spline shape functions [48].

This choice implies that the finite-dimensional spaces (37) and (38) for the
trial functions become

${{{}^{h/2}}\mathscr{V}_{\bar{\boldsymbol{u}}}}=\Bigl{\\{}{{}^{h/2}\boldsymbol{u}}\in\left[C^{m}\left({}^{h/2}\bar{\omega}\right)\right]^{n^{dim}}\bigm{|}\,{{}^{h/2}\boldsymbol{u}}\left(\boldsymbol{x}\right)=\mathrm{N}^{h/2}_{a}\left(\boldsymbol{x}\right)\mathbf{u}_{a}\\\
\forall
a\in\text{clos}\left({{}^{h/2}\mathcal{T}^{act}}\right),\,{{}^{h/2}\boldsymbol{u}}=\bar{\boldsymbol{u}}\
\hbox{on}\ {{}^{h/2}\bar{\gamma}^{\bar{u}}}\Bigr{\\}};$ (45)
${{{}^{h}}\mathscr{W}_{\bar{p}}}=\left\\{{}^{h}p^{(f)}\in
C^{m}\left({}^{h}\bar{\omega}\right)\bigm{|}\,{{}^{h}p^{(f)}}\left(\boldsymbol{x}\right)=\mathrm{N}^{h}_{a}\left(\boldsymbol{x}\right)\,\mathrm{p}_{a}\
\forall a\in\text{clos}\left({{}^{h}\mathcal{T}^{act}}\right)\right\\},$ (46)

and similar definitions follow for the spaces of test functions
${{}^{h/2}}\mathscr{V}_{\boldsymbol{0}}$ and
${{}^{h}}\mathscr{W}_{\boldsymbol{0}}$. Using QkSD-Qk-inspired meshes and
including the face ghost stabilisation Eq. (43), the weak form becomes: find
${{}^{h/2}\boldsymbol{u}}\in{{}^{h/2}\mathscr{V}_{\bar{\boldsymbol{u}}}}$ and
${{}^{h}p^{(f)}}\in{{{}^{h}}\mathscr{W}_{\bar{p}}}$ such that, for
$t\in[0,\tilde{t}]$,

$\int_{{{}^{MP}\omega}}\updelta\mathbf{w}_{a}\,\text{{grad}}\left(\mathrm{N}^{h/2}_{a}\right)\boldsymbol{:}\left(\boldsymbol{\sigma}^{\prime}-\mathrm{N}^{h}_{b}\,\mathrm{p}_{b}\,\boldsymbol{1}\right)\
dv\\\
-\int_{{{}^{MP}\omega}}\left(\updelta\mathbf{w}_{a}\,\mathrm{N}^{h/2}_{a}\right)\cdot\rho\boldsymbol{b}\
dv-\int_{\gamma^{\bar{t}}}\left(\updelta\mathbf{w}_{a}\,\mathrm{N}^{h/2}_{a}\right)\cdot\bar{\boldsymbol{t}}\
da\\\
+\gamma^{A}\,\frac{h}{2}\int_{{{}^{h/2}}\gamma^{F}}\updelta\mathbf{w}_{a}\left[\left[\left(\text{{grad}}\left(\mathrm{N}^{h/2}_{a}\right)\cdot\mathbf{n}\right)\right]\right]\cdot\left[\left[\left(\text{{grad}}\left(\mathrm{N}^{h/2}_{b}\right)\cdot\mathbf{n}\right)\right]\right]\mathbf{u}_{b}\
da=0,\\\ \quad\forall\updelta\mathbf{w}_{a}\in\
^{h/2}\mathscr{V}_{\boldsymbol{0}};$ (47)
$\int_{{{}^{MP}\omega}}\left(\updelta\upeta_{a}\,\mathrm{N}^{h}_{a}\right)\rho_{0}^{(f)}\frac{\dot{J}}{J}\
dv+\int_{{{}^{MP}\omega}}\frac{\upkappa}{g}\,\updelta\upeta_{a}\,\text{{grad}}\left(\mathrm{N}^{h}_{a}\right)\cdot\left(\text{{grad}}\left(\mathrm{N}^{h}_{b}\right)\mathrm{p}_{b}-\rho^{(f)}\boldsymbol{b}\right)dv\\\
-\int_{\gamma^{\bar{q}}}\left(\updelta\upeta_{a}\,\mathrm{N}^{h}_{a}\right)\,\bar{q}\
da+\gamma^{pen}\int_{\gamma^{\bar{p}}}\left(\updelta\upeta_{a}\,\mathrm{N}^{h}_{a}\right)\,\left(\mathrm{N}^{h}_{b}\,\mathrm{p}_{b}-\bar{p}\right)da\\\
+\,\gamma^{C}h\int_{{}^{h}\gamma^{F}}\updelta\upeta_{a}\left[\left[\left(\text{{grad}}\left(\mathrm{N}^{h}_{a}\right)\cdot\mathbf{n}\right)\right]\right]\left[\left[\left(\text{{grad}}\left(\mathrm{N}^{h}_{b}\right)\cdot\mathbf{n}\right)\right]\right]\mathrm{p}_{b}\
da=0,\quad\forall\updelta\upeta_{a}\in\ {{{}^{h}}\mathscr{W}_{0}},$ (48)

giving the initial condition
${{}^{h/2}\boldsymbol{u}}\left(t=0\right)=\boldsymbol{u}_{0}$. In the above
equations, $\gamma^{A}$ and $\gamma^{C}$ are the coefficients used to the
scale the effect of the ghost penalty terms and to match the physical units of
the equations. For this latter purpose, the units of $\gamma^{A}$ must be
similar to a stress measure, while the units of $\gamma^{C}$ should be the
same of the $\upkappa/g$.

Fig. 2(a) shows the application of the face ghost stabilisation to the finer
mesh of length $h/2$ for the displacement field (i.e., ghost penalty term
added to Eq. (47)), distinguishing between inactive and active elements,
boundary elements and boundary facets. Similarly, Fig. 2(b) illustrates the
same quantities for the coarser mesh of length $h$ employed for the pressure
field (i.e., ghost penalty term added to Eq. (48)). The combined result of the
different stabilised surfaces for the displacement and pressure fields are
represented in Fig. 2(c).

(a) Application of Algorithms 1 and 2 to the finer mesh for
${}^{h/2}\boldsymbol{u}$.

(b) Application of Algorithms 1 and 2 to the coarser mesh for ${}^{h}p^{(f)}$.

(c) Final result of facets where face ghost stabilisations are applied.

Figure 2: Algorithmic selection of boundary elements and facets (top row) for
finer (left) and coarser (right) meshes, and stabilised facets (bottom row).

For the sake of completeness, it must be acknowledged that a plethora of inf-
sup stabilised elements exist in the literature and covering the complete
literature on this subject is beyond the scope of this work. As for the FEM,
the reader can refer to, for instance, Leborgne [49] for a recent publication
discussing the inf-sup condition for different mixed formulations (and
possible remedies) or to the monograph of Boffi _et al._ [31] for a complete
discussion. The nature of QkSD-Qk elements makes its adaptation
straightforward to use in the MPM, especially if low-order shape functions are
used. On top of this, stable elements by design possess the advantage of not
adding terms to the primary equations. This is particularly advantageous if
the NR method is employed to solve the implicit primary equations, since less
terms can be considered in the Jacobian matrix. In the following section,
QkSD-Qk meshes are diversified for the sMPM and the GIMPM and their inf-sup
stability is investigated.

####  4.4.1 sMPM shape functions

Having selected a Cartesian mesh composed by quadrilateral (or hexahedral, if
$n^{dim}=3$) elements, the basis functions are constructed using the tensor
product of one-dimensional functions. In the case of sMPM for the generic mesh
of length $\tilde{h}$, first-order Lagrange polynomials are employed, i.e.,

$\displaystyle\mathrm{N}^{\tilde{h},1}\left(\xi\right)$
$\displaystyle=\begin{cases}1+\xi/\tilde{h},\qquad&\hbox{if}\quad-\tilde{h}<\xi\leq
0;\\\ 1-\xi/\tilde{h},\qquad&\hbox{if}\quad 0<\xi\leq\tilde{h},\end{cases}$
(49)

where $\xi={{}^{mp}x}-{{}^{\tilde{h}}x_{v}}$ denotes the one-dimensional local
coordinate computed as the difference between the material point’s coordinate
and the grid vertex coordinate. The use of these shape functions for a Q1SD-Q1
FEM $\boldsymbol{u}-p^{(f)}$ formulation has been tested in Dortdivanlioglu
_et al._ [30], employing the patch test proposed by Chapelle and Bathe [33]
for mixed FEM formulations. Employing this result and given the similarities
between the FEM and the MPM, it can be concluded that Q1SD-Q1 mesh-based
solutions for sMPM discretisation inherit inf-sup stability from its FEM
counterpart. This conclusion is also substantiated by Example 5.1.

####  4.4.2 GIMPM shape functions

The GIMPM shape functions in 1D are obtained via the convolution of a constant
function denoting the length of the particle $2l_{p}$ (named the
_characteristic function_) with the sMPM shape function given by Eq. (49). The
reader is referred to Bardenhagen _et al._ [36] or Charlton _et al._ [37] for
this convolution integral procedure. The one-dimensional GIMPM shape functions
defined on a mesh of length $\tilde{h}$ is given by

$\displaystyle\mathrm{S}^{\tilde{h}}\left(\xi\right)$
$\displaystyle=\begin{cases}\left(\tilde{h}+l_{p}+\xi\right)^{2}/\left(4\tilde{h}l_{p}\right),\qquad&\hbox{if}\quad-\tilde{h}-l_{p}<\xi\leq-\tilde{h}+l_{p};\\\
1+\xi/\tilde{h},\qquad&\hbox{if}\quad-\tilde{h}+l_{p}<\xi\leq-l_{p};\\\
1-\xi^{2}/\left(2\tilde{h}l_{p}\right)-l_{p}^{2}/\left(2\tilde{h}l_{p}\right),\qquad&\hbox{if}\quad-
l_{p}<\xi\leq l_{p};\\\ 1-\xi/\tilde{h},\qquad&\hbox{if}\quad
l_{p}<\xi\leq\tilde{h}-l_{p},\\\
\left(\tilde{h}+l_{p}-\xi\right)^{2}/\left(4\tilde{h}l_{p}\right),\qquad&\hbox{if}\quad\tilde{h}-l_{p}<\xi\leq\tilde{h}+l_{p},\end{cases}$
(50)

with $\xi={{}^{mp}x}-{{}^{\tilde{h}}x_{v}}$ as above.

Figure 3: Combination of GIMPM functions on overlapping meshes. Thick grey
lines illustrate the finer mesh, while white dashed lines represent the
coarser one.

Fig. 3 shows that three combinations of the polynomials are locally available
when using the GIMPM shape functions for QkSD-Qk meshes: Q2SD-Q2, Q1SD-Q1, and
Q2SD-Q1. As for the case of linear polynomials in Section 4.4.1, these
elements have been separately tested by Dortdivanlioglu _et al._ [30] for the
FEM using the patch test in Chapelle and Bathe [33]. When examined
independently, each of these elements has been shown to be inf-sup stable.
While it is not possible to apply this test directly to an MPM formulation
(patch tests cannot be adapted to the MPM), a thought experiment can be built
by employing shape functions Eq. (50) to construct a FEM discretisation and
using the results from Dortdivanlioglu _et al._ [30]. The test proposed by
Bathe and Chapelle [33] checks an eigenvalue-eigenvector problem using
matrices relative to discrete bilinear forms. This test is repeated for
different decreasing mesh sizes, and it is considered as passed when the
minimum eigenvalue does not decrease with finer grids. The same matrices for
the eigenvalue-eigenvector problem employing shape functions in Eq. (50) are a
linear combination of the above elements. This holds since each integration
point singularly contributes as Q2SD-Q2 or Q1SD-Q1 or Q2SD-Q1 to its element
and the assembling process sums each integration points’ contribution. Hence,
the boundedness of the minimum eigenvalue resulting from this FEM-like
discretisation using shape functions in Eq. (50) follows since it is linear
algebra (eigenvalue problem) applied to matrices of bilinear forms. This
rationale is confirmed by numerical Example 5.3, which adopts GIMPM shape
functions and exhibits stable pressure values.

###  4.5 Time discretisation

In addition to the spatial discretisation introduced above, a temporal
discretisation is required and this work considers the following Backward-
Euler time integration relative to the divergence of the velocity field

$\text{{div}}\cdot\boldsymbol{v}=\frac{\dot{J}}{J}=\left(\ln
J\right)^{\cdot}\approx\frac{\left(\ln J_{n+1}-\ln J_{n}\right)}{\Delta t},$
(51)

where the subscript $n+1$ denotes the current time-step and $n$ the previous
one.

The literature (see, for instance, Sun _et al._ [50] or Zhao and Choo [51])
recognises the role played by the above time-discretisation to avoid negative
values of the Jacobian. However, the current formulation considering an
incompressible solid phase imposes a more severe constraint on the Jacobian
than $J>0$. This can quickly verified if inequalities (3) relative to the
Eulerian porosity are substituted in solid mass conservation Eq. (8). The
reader can refer to Pretti _et al._ [28] for further details and consequences
of the above inequalities. Despite this consideration, the above formula
possesses another desirable feature for the MPM, discussed below in Section
4.6 and it is thus employed.

###  4.6 Mapping processes and MPs’ domain update

As mentioned in Section 4.1, two mapping processes take place in an MPM step.
The Point-to-Grid (P2G) mapping (phase $(\mathrm{ii})$ in Figure 1)
initialises the nodal grid unknowns, transferring the information from the
MPs. This mapping process is required to conserve physical properties of
interest relative to the mapped information. In the context of the dynamics of
a solid body, for instance, the velocity mapped from the MPs to the grid is
expected to conserve, as much as possible, momenta (linear and angular) and
kinetic energy to avoid numerical dissipation (see Love and Sulsky [52] or
Pretti _et al._ [53]). A similar rigorous procedure should be expected for a
poro-mechanical problem. However, mapping the velocity can be avoided for a
quasi-static problem using Eq. (51), since the Jacobian $J$ is not required on
the grid for computational purposes. As for the other primary variable, the
fluid pressure relates to a form of energy stored in the fluid body only if
the considered material is compressible, as pointed out by Miehe _et al._
[29]. Owing to Assumption (6), this work considers an incompressible fluid
constituent, which is not required to be mapped as the pressure acts as a
Lagrange multiplier, which does not entail any form of energy. Overall, it
emerges that G2P mapping is not required in this context, i.e., quasi-static
simulations considering an incompressible fluid phase. However, this procedure
can take place for estimating better predictors for the NR scheme, but these
estimates are free from numerical dissipation .

Once the grid solution is achieved, this must be mapped to the MPs before grid
disposal (phase $(\mathrm{v})$ in Figure 1). Unlike the P2G mapping, the
Grid-2-Point (G2P) mapping must always take place, as it is the MPs’ role to
store information necessary to run the analysis through different steps. Under
the same assumptions considered in the above paragraph for the P2G mapping
(i.e., quasi-static simulations involving an incompressible fluid), the G2P
mapping is not required to conserve physical quantities of interest for the
same reasons outlined above. Hence, the shape functions can be
straightforwardly used to map the values of the unknowns from the grid nodes
to the MPs, i.e.,

$\displaystyle{{}^{mp}\boldsymbol{u}_{n+1}}={{}^{mp}\boldsymbol{u}_{n}}+\mathrm{N}^{h/2}_{a}\left({}^{mp}\boldsymbol{x}_{n+1}\right)\,\Delta\mathbf{u}_{a};$
(52)
$\displaystyle{{}^{mp}\boldsymbol{p}_{n+1}}=\mathrm{N}^{h}_{a}\left({}^{mp}\boldsymbol{x}_{n+1}\right)\,\left(\mathbf{p}_{n+1}\right)_{a},$
(53)

where $\Delta\mathbf{u}_{a}$ is the difference in the time-step of the grid
displacements, i.e.,
$\Delta\left(\bullet\right)\coloneqq\left(\bullet\right)_{n+1}-\left(\bullet\right)_{n}$.
The above updates exploit the considerations drawn so far: on the one hand,
the initial displacement on the grid $\left(\mathbf{u}_{a}\right)_{n}$ is not
reconstructed with the P2G mapping (the unknown of the system of primary
equations can directly be $\Delta\mathbf{u}_{a}$); on the other, the pressure
field can be re-written both at the MP and grid level because it is not
related to an energy measure.

While the update of the mixed volume for the sMPM follows the standard
procedure for volume update using the Jacobian, the length of the
characteristic function for the GIMPM requires an update which can differ for
the different Cartesian directions. This work follows the _corner update
procedure_ proposed by Coombs _et al._ [22]. The reader is referred to this
reference for further details on its implementation, and the exhibited
advantages over the other techniques.

##  5 Numerical examples

The model detailed in Section 2, formulated in a weak sense in 3, and
discretised in 4 is below investigated via three different examples, each
exploring a particular feature of the proposed formulation. The implementation
of the outlined model has been carried out in an extended version of AMPLE
[54].

###  5.1 Terzaghi mono-dimensional consolidation

##### Example scope

The numerical investigation provided by the Terzaghi mono-dimensional
consolidation (see, for instance, [55]) delivers a two-fold goal. A comparison
considering two meshes is made between the QkSD-Qk meshes (specifically using
sMPM, i.e., Q1SD-Q1) and the Polynomial Pressure Projection (PPP) with a
single mesh. The PPP is one of the most widely adopted stabilisation
techniques in the literature for $\boldsymbol{u}-p^{(f)}$ (see White and Borja
[56] for the FEM, and Zhao and Choo [51] for its adaptation to the MPM). The
results on the two meshes are compared near zero consolidation time versus an
established critical time for $\boldsymbol{u}-p^{(f)}$ formulations (see
Vermeer [57]). This time sets that the coarser simulation shows oscillating
pressure, while the finer analysis does not. While the simulations with
Q1SD-Q1 meshes agree with these predictions, it is shown that the analyses
employing the PPP exhibits pressure instability for both meshes. Furthermore,
in the case of the coarser unstable simulation, it is investigated how this
instability dissipates over time for the simulation using Q1SD-Q1 meshes.

The analytical values of pressure are well-known in the small-strain regime
and are given by (see, for instance, [55])

$P(Z,T)=\sum_{m=0}^{\infty}\frac{2}{M}\,\sin\left(MZ\right)\exp\left(-M^{2}\,T\right),\\\
\hbox{with}\ M=\frac{\pi}{2}\left(2m+1\right)\quad\text{and}\
m=0,1,\dots\in\mathbb{N}_{0},$ (54)

where the adimensional quantities in the above equation are defined as follows

$P\coloneqq\frac{p^{(f)}}{w};\qquad Z\coloneqq\frac{z}{H};\qquad
T\coloneqq\frac{c_{v}}{H^{2}}\,t,$ (55)

with $w$ being the magnitude of the overburden and $H$ the height of the
column. The relationship between the coefficient of consolidation $c_{v}$ and
other hydro-mechanical parameters is given by

$c_{v}=\frac{\left(\frac{K}{n_{0}}+\frac{4}{3}G\right)}{\rho^{(f)}\,g}\,\upkappa,$
(56)

where $\frac{K}{n_{0}}$ is the initial tangent modulus of the improved Hencky
material described by Eq. (15).

Table 1: Summary of the parameters considered in the analysis of the Terzaghi
mono-dimensional consolidation. Parameter Settings
---
Material Parameters | $\frac{K}{n_{0}},\,G$ |  $1\cdot 10^{6}$ Pa, $6\ \cdot 10^{5}$ Pa
$\upkappa_{0}$ |  $10^{-7}$ m s-1
$\rho^{(f)}$ |  1000 kg/m3
Geometry and loading | $H$ |  $1$ m
$w$ |  $100$ kPa

Figure 4: Illustration of the Terzaghi mono-dimensional problem. Permeable
surfaces are designed by the dashed line.

##### Setup

The hydro-mechanical parameters are the same employed by Zhao and Choo [51]
and are reported in Table 1. Owing to the assumed parameters, it is expected
that the numerical values, computed in the finite strain context, can
reproduce the analytical results set in the small strain theory. Two MPs per
direction are initially positioned for each element of the finer mesh. Owing
to this initial discretisation and given that small displacements are
expected, the small-cut issue detailed in Section 4.2 does not occur in this
example. Thus, the ghost stabilisation included in Eqs. (47) and (48) are not
considered for this problem.

As mentioned above, two different vertical discretisations have been
considered, and the initial selected time-step is $t_{0}=0.1$ s. These
discretisations have been designed to comply with (in the case of
$n_{y}^{els}=320$) and violate (in the case of $n_{y}^{els}=160$) at the first
time-step the time value prescribed by Vermeer and Verruijt [57] for implicit
$\boldsymbol{u}-p^{(f)}$ formulations, this being

$t\geq\frac{\Delta h_{y}^{2}}{6\,c_{v}}\begin{cases}\approx
0.09\,\text{s}\qquad\text{if}\ n^{els}_{y}=320;\\\ \approx
0.36\,\text{s}\qquad\text{if}\ n^{els}_{y}=160.\end{cases}$ (57)

To avoid unnecessary computational time and still provide a good description
of the consolidation phenomenon, the time-step partition observes the
geometric series to complete the consolidation process (i.e., $T=1$)

$\sum_{p=0}^{n-1}c^{p}t_{0}\approx\frac{H^{2}}{c_{v}},$ (58)

where the selected common ratio $c=1.01673$ has been selected to run the
simulation in $n=550$ time-steps.

The analyses have been run using sMPM shape functions (Q1SD-Q1 meshes).
However, given that small displacements are expected, no appreciable
difference could be noticed if the GIMPM shape function should have been
employed.

The parameter for the PPP vector is $\frac{1}{2G}$, which is widely adopted in
the literature (see Zhao and Choo [51] and references therein).

##### Results discussion

(a) Excess pore pressure isochrones comparison at the time $T\approx 1.8\cdot
10^{-6}$ s for the case with $n^{els}_{y}=160$.

(b) Excess pore pressure isochrones comparison at the time $T\approx 1.8\cdot
10^{-6}$ s for the case with $n^{els}_{y}=320$.

Figure 5: Graphical comparison between analytical formula and numerical
solutions obtained with the PPP stabilisation and Q1SD-Q1 meshes.

From a comparison between Figs. 5(a) and 5(b), it is clear that the
formulation proposed in this paper (green line) agrees with the prediction of
unstable (Fig. 5(a)) and stable (Fig.5(b)) time provided by Eq. (57).
Conversely, the PPP (red line) exhibits unstable behaviour even for the
predicted stable discretisation with $n_{y}^{els}=320$. Moreover, it can be
seen that the instability caused by the PPP propagates through more vertical
elements than its Q1SD-Q1 counterpart for both situations. It must be observed
that these results comply with the literature: Preisig and Prevost [58] stated
that the PPP is unable to remove pressure oscillation near the draining
boundary, and Monforte _et al._ [59], adopting a
$\boldsymbol{u}-\boldsymbol{w}-p^{(f)}$ formulation, showed pressure
oscillating behaviour for a high-frequency wave propagation problem. These
studies confirm that the PPP was not designed to stabilise poro-mechanical
problems in which the pressure field belongs to $H^{1}(\omega)$: the original
application of the PPP (see Dohrmann and Bochev [60]) was Stokes problems,
which requires the pressure field to be $H^{0}(\omega)$. The peaks shown in
Figs. 5(a) and 5(b) (and confirmed by the literature) correspond with sharp
pressure gradients (i.e., with the phenomenon progressively passing from
undrained to drained), which the PPP was not designed to stabilise.

Figure 6: Excess pore pressure isochrones for the case with $n^{els}_{y}=160$.

Fig. 6 (particularly its magnification) investigates the persistence of the
spurious peak for the Q1SD-Q1 formulation in the unstable case (i.e.,
$n^{els}_{y}=160$). The numerical results agree with Eq. (57), which
prescribes a stable time bigger than $\approx 0.36$ s. For the adimensional
time $T\approx 1.33\cdot 10^{-5}$ (corresponding to $\approx 0.74$ s) the peak
is practically extinguished. Its pressure value is $2.14\%$ bigger than its
analytical one, which is significantly lower than the one exhibited in Fig.
5(a) for $T\approx 1.8\cdot 10^{-6}$ (corresponding to $\approx 0.1$ s), this
being $17.5\%$ bigger compared to the analytical value. For the $T\approx
1.05\cdot 10^{-4}$ (corresponding to $\approx 5.83$ s), there is no
appreciable difference between the numerical and the analytical solution. This
overlap between the numerical investigation and the analytical solution for
the exceed pore pressure isochrones continues for the rest of the simulation,
until consolidation has entirely taken place.

###  5.2 Investigation of the face ghost penalty

Figure 7: Geometry and initial setup for the considered ellipse. Each initial
setup consists in seeding the MPs discretising the ellipse in a different
position, described the coordinate $a$, while keeping the meshes (finer in
grey, coarser dashed in white) fixed. New MP-setups are generated for
different values of $a$. Drained conditions are applied at the top of the mesh
(dashed black line).

##### Example scope

While the previous example focused on the stabilising effect given by Q1SD-Q1
elements, this numerical exercise (illustrated in Figure 7) is particularly
aimed at investigating the beneficial role of the ghost penalties when applied
to the matrices $\mathbf{A}$ and $\mathbf{C}$ appearing in the Jacobian
matrix. If adopting (as is the case in this example) QkSD-Qk meshes, the
drained case requires coercivity of both these bilinear forms associated with
the above matrices to guarantee solvability.

The effectiveness of the face ghost stabilisation applied to Eqs. (47) and
(48) is investigated for different values of the parameters $\gamma^{A}$ and
$\gamma^{C}$. The spanned values of these parameters have been chosen to be
not too small, resulting in insignificant effects, nor too high, leading to
potential locking phenomena (see Badia et al. [61]). This example bears
similarities with other numerical tests proposed by Coombs [40] in the MPM
context and by Kothari and Krause [62] and Sticko _et al._ [63] for the
unfitted FEM.

Parameter Settings
---
Material Parameters | $\frac{K}{n_{0}},\,G$ |  $1$ Pa, $1$ Pa
$\upkappa_{0}$ |  $1$ m s-1
$\rho^{(f)}$ |  $1$ kg/m3
Geometry | $H_{x},\,H_{y}$ |  $8$ m, $3$ m
$r_{x},\,r_{y}$ |  $2$ m, $2.15$ m
Grid | $\Delta h_{x},\,\Delta h_{y}$ |  $1$ m, $1$ m
Table 2: Summary of the parameters considered in the investigation of the
ghost penalty example.

##### Setup

The half of an ellipse illustrated in Fig. 7 is made of an ideal porous
material, whose parameters are listed in Table 2. Since the scope of this
investigation lies in the examination of the small-cut issue (and its remedy),
the geometrical aspect, i.e., the overlap between the physical domain and
shape functions’ stencil, is examined. This justifies two aspects: on the one
hand, the assumption of ideal, non-descriptive, hydro-mechanical parameters;
on the other, the use of a geometry such an ellipse. This form combines a
periodic behaviour along the coordinate $a$ (see Fig. 7) and aggravates the
small-cut problem at the bottom limit of the ellipse.

Three MPs per direction are initially equally distributed per each element of
the finer mesh. The MPs lying outside the analytical shape described by the
half of the ellipse are then removed. This setup is repeated for different
values of the coordinate $a$, which is varied progressively for each
simulation. For each of these discretisations, the submatrices $\mathbf{A}$
and $\mathbf{C}$ in the Jacobian are assembled, so that the effect of the
ghost stabilisation can be assessed for different physical MP-based domain-
mesh interactions.

It must be noted that, for the sMPM, it is not possible to assign an element
to a MP when this lie precisely on an element boundary. Since this occurs for
the finer mesh when $a=\frac{\Delta h_{x}}{4}*p$, with
$p=1,2,\dots\in\mathbb{N}$, these sampling location are avoided for the sMPM.
This choice explains the rupture of periodic behaviour which is particularly
evident in Figure 8(a). Nonetheless, it must be underlined that this situation
is due to the specificity of the setup and is highly unlikely to occur in
standard simulations.

(a) Condition numbers for the $\mathbf{A}$ sub-matrix for translating domains,
sMPM.

(b) Condition numbers for the $\mathbf{C}$ sub-matrix for translating domains,
sMPM.

(c) Condition numbers for the $\mathbf{A}$ sub-matrix for translating domains,
GIMPM.

(d) Condition numbers for the $\mathbf{C}$ sub-matrix for translating domains,
GIMPM.

Figure 8: Effects of the applied ghost penalty parameters on the condition
number of submatrices $\mathbf{A}$ (left column) and $\mathbf{C}$ (right
column) in the case of sMPM (top row) and GIMPM. (bottom row).

##### Results discussion

As it can be appreciated from Fig. 8, there is considerable difference in the
condition numbers (denoted by $\upkappa(\bullet)$, with $\bullet$ being the
considered submatrix) for unstabilised sMPM and GIMPM, with the latter showing
on average three orders of magnitude higher condition number for the same
setup. This behaviour is due to the shape functions stencil, which is more
extended for the GIMPM.

Fig. 8(a) highlights how the lowest value of ghost penalty ($\gamma^{A}=E\cdot
10^{-4}$) does not contribute significantly to the submatrix $\mathbf{A}$
condition number in the case of sMPM. Diversely, in the case of the GIMPM in
Fig. 8(c), the effect of the ghost penalty for $\gamma^{A}=E\cdot 10^{-4}$
reduces the peak values by approximately an order of magnitude, even though
the general pattern follows similarly that of the unstabilised version. These
peaks are more levelled for the case of $\gamma^{A}=E\cdot 10^{-2}$ and
entirely smoothed for the case of $\gamma^{A}=E$, which exhibits condition
numbers approximately 4.5 order of magnitude lower than the peaks of the
unstabilised version. Back to Fig. 8(a) for the sMPM, it can be seen how
$\gamma^{A}=E$ leads to overall smooth behaviour, with the condition number
being approximately 1.5 orders of magnitude below the unstabilised peaks. The
value $\gamma^{A}=E\cdot 10^{-2}$ contributes less significantly to
stabilising the sMPM than its GIMPM counterpart. Similar trends can be seen
for the condition number of the submatrix $\mathbf{C}$ in the case of sMPM
(Fig. 8(b)) and GIMPM (Fig. 8(d)). This behaviour is interesting, as the
submatrix $\mathbf{C}$ is entirely assembled using the coarser grid of length
$h$. Hence, it highlights how coarser discretisation, which should be less
prone to the small-cut issue when using the same numbers of MPs, is not free
from this problem and strengthens the need for a stabilisation such as the one
included in this work.

Owing to the previously mentioned lack of rigorous analysis for the MPM, the
condition numbers of the submatrices $\mathbf{A}$ and $\mathbf{C}$, even if
they benefit from the ghost stabilisation, are still dependent to some degree
on the intersection between the shape functions’ stencil and the MP-based
physical domain. This aspect diversifies the MPM from the unfitted FEM, where,
once stabilised, the domain’s cuts do not play a role in the condition number
(as proved in Burman [39] and numerically demonstrated in Kothari and Krause
[62]). This difference amplifies the importance of carefully choosing the
correct value of the stabilisation parameter, as some selections, even though
effective, might not give the desired limitation on the condition number.

###  5.3 Flexible strip foundation

Parameter Settings
---
Material Parameters | $K,\,G$ |  $1.062\cdot 10^{8}$ Pa, $9.8\cdot 10^{7}$ Pa
$\alpha,\,\gamma,\,p_{c}$ |  $0.6,\,0.9,\,7.5\cdot 10^{7}$ Pa
$\upkappa_{0}$ |  $2.07\cdot 10^{-3}$ m s-1
$n_{0}$ | 0.5
$\rho^{(f)}$ |  $1000$ kg/m3
Geometry and loading | $L,\,H,\,l$ |  $20$ m, $16$ m, $8$ m
$w$ |  $2.5\cdot 10^{8}$ Pa
Time and time-step | T |  $10$ s
$\Delta t$ |  $0.5$ s
Grid and MPs | $\Delta h_{x},\,\Delta h_{y}$ |  $0.25$ m, $0.25$ m
MPs per element (finer mesh) | 4
Table 3: Summary of the parameters considered in the analysis of the flexible
strip foundation. Figure 9: Illustration of the flexible strip foundation
problem. Permeable surfaces are designed by the dashed line.

##### Example scope

The scope of this numerical test is to support the rationale provided in
Section 4.4.2 regarding the stability of the different combinations of
polynomial functions employed by GIMPM formulation with QkSD-Qk meshes. To
demonstrate that the pressure field is free from oscillations, a flexible
footing loading has been applied at the top of a fully saturated porous
material in plane-strain conditions. Elastic and the elasto-plastic cases have
been considered. Since the application of the loading causes the considered
body to displace substantially, it is expected that all the combinations of
polynomials that the GIMPM basis functions exhibit will be met through the
analysis, even if not uniformly or simultaneously.

##### Setup

Hydro-mechanical parameters are reported in Table 3, and the considered
elastic material shares similarities with the one presented in Armero [64] in
a very similar example. The modified Cam clay $\alpha-\gamma$ yield surface
proposed by Collins and Hilder [65] has been considered for the elasto-plastic
case (see, for instance, [66] for implementation details).

To apply the non-conforming pressure BCs, a penalty value of
$\gamma^{pen}=5\cdot 10^{6}/\upkappa$ was employed, where the hydraulic
conductivity appears at the denominator for a matter of consistency with the
physical units. The load has been linearly increased from zero to the value
shown in Table 3. The elastic case can complete the simulation in 20 time-
step, while the plastic case reaches the limit load in correspondence of the
$4^{th}$ time-step, where it is stopped.

(a) Contours of the fluid pressure field, $1^{st}$ time-step.

(b) Contours of the fluid pressure field, $5^{th}$ time-step.

(c) Contours of the fluid pressure field, $10^{th}$ time-step.

(d) Contours of the fluid pressure field, $15^{th}$ time-step.

(e) Contours of the fluid pressure field, $20^{th}$ time-step.

Figure 10: Contours of the pressure field across the simulation for the
flexible footing loading, elastic case.

Given the low number of MPs per element and the substantial displacements
expected in this simulation, the ghost penalisation have been included. Based
on the observations made in Example 5.2, the selected parameter for the
penalty value are $\gamma^{A}=E\cdot 10^{-1}$ and
$\gamma^{C}=\frac{\upkappa}{g}\cdot 10^{-1}$.

##### Results discussion

As can be seen from the contours shown throughout the analysis in Fig. 10 for
the elastic case, the fluid pressure does not present oscillations typical of
violating the inf-sup condition. Even though showing spatially different
distributions of pressure, these values are equally smooth for the final step
of the elasto-plastic case shown in Fig. 11(b). As Fig. 11(a) shows, plastic
strains arise in correspondence with the footing angle, forming the expected
wedge shape. This shape is also visible for the pressure field in Fig. 11(b),
where the GIMPM QkSD-Qk meshes handle this horizontal gradient without
interruptions or oscillations. Due to the continuity demonstrated for the
different fields and situations in Figs. 10 and 11, the chosen ghost
parameters have proved suitable for stabilising the submatrices’ conditioning
number and achieving the continuous expected results.

Overall, the exhibited smooth values of pressure substantiate the explanation
provided in Section 4.4.2 for the GIMPM shape functions, showing that these
described QkSD-Qk stable meshes deliver oscillating-free solutions.

(a) Contours of the plastic volumetric part of the logarithmic strain, last
time-step.

(b) Contours of the fluid pressure field, last time-step.

Figure 11: Contours of different fields at the last time-step for the flexible
footing loading, plastic case.

##  6 Conclusions

This manuscript has highlighted two sources of instability that can arise for
the MPM in the case of mixed formulations and proposed a new approach that
overcomes both issues for coupled (solid-fluid) problems. The intrinsic nature
of the MPM as an unfitted method and its resulting small-cut issue can lead to
ill-conditioned matrices. For mixed formulations, this small-cut issue affects
both the submatrices appearing on the main diagonal of the Jacobian (i.e.,
$\mathbf{A}$ and $\mathbf{C}$), especially in drained conditions. The
treatment of this instability has been the use of the face ghost penalty
method on both the displacement and fluid pressure meshes, which, even though
it does not guarantee coercivity for the MPM, limits the condition number of
these submatrices, thus permitting the inversion of the Jacobian matrix.

Furthermore, in nearly undrained conditions, the choices of the spaces of test
and trial functions relative to displacement and pressure can violate the inf-
sup condition. This work has adapted the use of QkSD-Qk elements to the MPM.
The resulting overlapping meshes, while maintaining low-order of shape
functions (as is the case with the sMPM and the GIMPM), are stable by design.
This is in contrast to the widely adopted PPP formulation, which is unable to
stabilise poro-mechanical problems, especially near drained boundaries, due to
the $H^{1}(\omega)$ nature of the pressure field [58].

The key contribution offered by this paper is a stable, implicit MPM
formulation for large deformation coupled problems including elasto-plastic
material behaviour. While a rigorous analysis or even the patch test are not
available for the MPM, the rationales provided throughout the paper are
substantiated via the similarities between the MPM and the unfitted FEM and
tested with the considered numerical examples.

## Acknowledgements

The first author was supported by the Faculty of Science at Durham University.
The second, fourth and fifth were supported by Engineering and Physical
Sciences Research Council, Grant Number EP/W000970/1. The third author was
supported by Physical Sciences Research Council, Grant Number EP/T518001/1.
All data created during this research are openly available at
collections.durham.ac.uk/ (specific DOI to be confirmed if/when the paper is
accepted). The research presented in this article has also benefited from
discussions with, and feedback from Ted O’Hare, Bradley Sims, Sam Sutcliffe,
and Dr. Mao Ouyang. For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising.

## References

  * [1] D. Sulsky, Z. Chen, H. Schreyer, A particle method for history-dependent materials, Computer Methods in Applied Mechanics and Engineering 118 (1-2) (1994) 179–196.
  * [2] D. Sulsky, S.-J. Zhou, H. L. Schreyer, Application of a particle-in-cell method to solid mechanics, Computer Physics Communications 87 (1-2) (1995) 236–252.
  * [3] J. Gaume, A. van Herwijnen, T. Gast, J. Teran, C. Jiang, Investigating the release and flow of snow avalanches at the slope-scale using a unified model based on the material point method, Cold Regions Science and Technology 168 (2019) 102847.
  * [4] J. Gaume, A. M. Puzrin, Mechanisms of slab avalanche release and impact in the dyatlov pass incident in 1959, Communications Earth & Environment 2 (1) (2021) 10.
  * [5] D. Sulsky, H. Schreyer, K. Peterson, R. Kwok, M. Coon, Using the material-point method to model sea ice dynamics, Journal of Geophysical Research: Oceans 112 (C2) (2007).
  * [6] A. Huth, R. Duddu, B. Smith, A generalized interpolation material point method for shallow ice shelves. 1: Shallow shelf approximation and ice thickness evolution, Journal of Advances in Modeling Earth Systems 13 (8) (2021) e2020MS002277.
  * [7] M. Xie, P. Navas, S. Lopez-Querol, A stabilised semi-implicit double-point material point method for soil-water coupled problems, arXiv preprint arXiv:2401.11951 (2024).
  * [8] N. D. Gavin, G. Pretti, W. M. Coombs, J. C. Brigham, C. E. Augarde, On the implementation of a material point-based arc-length method, International Journal for Numerical Methods in Engineering (2024) e7438.
  * [9] J. Davy, P. Lloyd, J. H. Chandler, P. Valdastri, A framework for simulation of magnetic soft robots using the material point method, IEEE Robotics and Automation Letters (2023).
  * [10] M. Li, Y. Lei, D. Gao, Y. Hu, X. Zhang, A novel material point method (mpm) based needle-tissue interaction model, Computer Methods in Biomechanics and Biomedical Engineering 24 (12) (2021) 1393–1407.
  * [11] A. Stomakhin, C. Schroeder, L. Chai, J. Teran, A. Selle, A material point method for snow simulation, ACM Transactions on Graphics (TOG) 32 (4) (2013) 1–10.
  * [12] J. Chen, V. Kala, A. Marquez-Razon, E. Gueidon, D. A. Hyde, J. Teran, A momentum-conserving implicit material point method for surface tension with contact angles and spatial gradients, ACM Transactions on Graphics (TOG) 40 (4) (2021) 1–16.
  * [13] O. Coussy, Poromechanics, John Wiley & Sons, 2004.
  * [14] X. Zheng, S. Wang, F. Yang, J. Yang, Material point method simulation of hydro-mechanical behaviour in two-phase porous geomaterials: A state-of-the-art review, Journal of Rock Mechanics and Geotechnical Engineering (2023).
  * [15] K. Soga, E. Alonso, A. Yerro, K. Kumar, S. Bandara, Trends in large-deformation analysis of landslide mass movements with particular emphasis on the material point method, Géotechnique 66 (3) (2016) 248–273.
  * [16] H. A. Nøst, Iteratively coupled implicit dynamic MPM-FVM, Ph.D. thesis, Norwegian University of Science and Technology (2019).
  * [17] X. Zheng, Stabilised material point method for fluid-saturated geomaterials, Ph.D. thesis, Delft University of Technology (2022).
  * [18] G. Pretti, Continuum mechanics and implicit material point method to underpin the modelling of drag anchors for cable risk assessment, Ph.D. thesis, Durham University (2024).
  * [19] E. Burman, P. Hansbo, Fictitious domain methods using cut elements: Iii. a stabilized nitsche method for stokes’ problem, ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 48 (3) (2014) 859–874.
  * [20] V. P. Nguyen, A. de Vaucorbeil, S. Bordas, The Material Point Method: Theory, Implementations and Applications, Springer Nature, 2023.
  * [21] J. E. Guilkey, J. A. Weiss, Implicit time integration for the material point method: Quantitative and algorithmic comparisons with the finite element method, International Journal for Numerical Methods in Engineering 57 (9) (2003) 1323–1338.
  * [22] W. M. Coombs, C. E. Augarde, A. J. Brennan, M. J. Brown, T. J. Charlton, J. A. Knappett, Y. G. Motlagh, L. Wang, On lagrangian mechanics and the implicit material point method for large deformation elasto-plasticity, Computer Methods in Applied Mechanics and Engineering 358 (2020) 112622.
  * [23] E. Kröner, Allgemeine kontinuumstheorie der versetzungen und eigenspannungen, Archive for Rational Mechanics and Analysis 4 (1) (1959) 273–334.
  * [24] E. H. Lee, Elastic-plastic deformation at finite strains, Journal of Applied Mechanics 36 (1) (1969) 1–6.
  * [25] J. Mandel, Plasticité Classique et Viscoplasticité, Vol. 97 of CISM Courses and Lectures, Springer, 1971.
  * [26] J. F. Thorpe, On the momentum theorem for a continuous system of variable mass, American Journal of Physics 30 (9) (1962) 637–640.
  * [27] R. I. Borja, E. Alarcón, A mathematical framework for finite strain elastoplastic consolidation part 1: Balance laws, variational formulation, and linearization, Computer Methods in Applied Mechanics and Engineering 122 (1-2) (1995) 145–171.
  * [28] G. Pretti, W. M. Coombs, C. E. Augarde, M. M. Puigvert, J. A. R. Gutiérrez, Preserving non-negative porosity values in a bi-phase elasto-plastic material under terzaghi’s effective stress principle, Mechanics of Materials (2024) 104958.
  * [29] C. Miehe, S. Mauthe, S. Teichtmeister, Minimization principles for the coupled problem of darcy–biot-type fluid transport in porous media linked to phase field modeling of fracture, Journal of the Mechanics and Physics of Solids 82 (2015) 186–217.
  * [30] B. Dortdivanlioglu, A. Krischok, L. Beirão da Veiga, C. Linder, Mixed isogeometric analysis of strongly coupled diffusion in porous materials, International Journal for Numerical Methods in Engineering 114 (1) (2018) 28–46.
  * [31] D. Boffi, F. Brezzi, M. Fortin, et al., Mixed finite element methods and applications, Vol. 44, Springer, 2013.
  * [32] P. Mira, M. Pastor, T. Li, X. Liu, A new stabilized enhanced strain element with equal order of interpolation for soil consolidation problems, Computer methods in applied mechanics and engineering 192 (37-38) (2003) 4257–4277.
  * [33] D. Chapelle, K.-J. Bathe, The inf-sup test, Computers & structures 47 (4-5) (1993) 537–545.
  * [34] B. Nedjar, Formulation of a nonlinear porosity law for fully saturated porous media at finite strains, Journal of the Mechanics and Physics of Solids 61 (2) (2013) 537–556.
  * [35] B. Nedjar, On finite strain poroplasticity with reversible and irreversible porosity laws. formulation and computational aspects, Mechanics of Materials 68 (2014) 237–252.
  * [36] S. G. Bardenhagen, E. M. Kober, The generalized interpolation material point method, Computer Modeling in Engineering and Sciences 5 (6) (2004) 477–496.
  * [37] T. Charlton, W. Coombs, C. Augarde, iGIMP: An implicit generalised interpolation material point method for large deformations, Computers & Structures 190 (2017) 108–125.
  * [38] R. E. Bird, G. Pretti, W. M. Coombs, C. E. Augarde, An implicit material point-to-rigid body contact approach for large deformation soil-structure interaction, Submitted to Computers and Geotechnics (2024).
  * [39] E. Burman, Ghost penalty, Comptes Rendus. Mathématique 348 (21-22) (2010) 1217–1220.
  * [40] W. M. Coombs, Ghost stabilisation of the material point method for stable quasi-static and dynamic analysis of large deformation problems, International Journal for Numerical Methods in Engineering 124 (21) (2023) 4841–4875.
  * [41] E. Burman, P. Hansbo, M. G. Larson, On the design of locking free ghost penalty stabilization and the relation to cutfem with discrete extension, arXiv preprint arXiv:2205.01340 (2022).
  * [42] Z. Liu, Y. Zhang, Y. Jiang, H. Yang, Y. Yang, Unfitted finite element method for fully coupled poroelasticity with stabilization, Computer Methods in Applied Mechanics and Engineering 397 (2022) 115132.
  * [43] Y. Yamaguchi, S. Moriguchi, K. Terada, Extended b-spline-based implicit material point method, International Journal for Numerical Methods in Engineering 122 (7) (2021) 1746–1769.
  * [44] K. Höllig, U. Reif, J. Wipper, Weighted extended b-spline approximation of dirichlet problems, SIAM Journal on Numerical Analysis 39 (2) (2001) 442–462.
  * [45] K. Höllig, Finite element methods with B-splines, SIAM, 2003.
  * [46] S. Badia, F. Verdugo, A. F. Martín, The aggregated unfitted finite element method for elliptic problems, Computer Methods in Applied Mechanics and Engineering 336 (2018) 533–553.
  * [47] M. Bercovier, O. Pironneau, Error estimates for finite element method solution of the stokes problem in the primitive variables, Numerische Mathematik 33 (1979) 211–224.
  * [48] A. A. Madadi, B. Dortdivanlioglu, A subdivision-stabilized b-spline mixed material point method, Computer Methods in Applied Mechanics and Engineering 418 (2024) 116567.
  * [49] G. Leborgne, Inf-sup condition and locking: Understanding and circumventing. stokes, laplacian, bi-laplacian, kirchhoff–love and mindlin–reissner locking type, boundary conditions, arXiv preprint arXiv:2301.04373 (2023).
  * [50] W. Sun, J. T. Ostien, A. G. Salinger, A stabilized assumed deformation gradient finite element formulation for strongly coupled poromechanical simulations at finite strain, International Journal for Numerical and Analytical Methods in Geomechanics 37 (16) (2013) 2755–2788.
  * [51] Y. Zhao, J. Choo, Stabilized material point methods for coupled large deformation and fluid flow in porous materials, Computer Methods in Applied Mechanics and Engineering 362 (2020) 112742.
  * [52] E. Love, D. L. Sulsky, An unconditionally stable, energy–momentum consistent implementation of the material-point method, Computer methods in applied mechanics and engineering 195 (33-36) (2006) 3903–3925.
  * [53] G. Pretti, W. M. Coombs, C. E. Augarde, B. Sims, M. M. Puigvert, J. A. R. Gutiérrez, A conservation law consistent updated lagrangian material point method for dynamic analysis, Journal of Computational Physics 485 (2023) 112075.
  * [54] W. M. Coombs, C. E. Augarde, AMPLE: a material point learning environment, Advances in Engineering Software 139 (2020) 102748.
  * [55] K. Terzaghi, Theoretical soil mechanics, John Wiley & sons, New York, 1943.
  * [56] J. A. White, R. I. Borja, Stabilized low-order finite elements for coupled solid-deformation/fluid-diffusion and their application to fault zone transients, Computer Methods in Applied Mechanics and Engineering 197 (49-50) (2008) 4353–4366.
  * [57] P. Vermeer, A. Verruijt, An accuracy condition for consolidation by finite elements, International Journal for numerical and analytical methods in geomechanics 5 (1) (1981) 1–14.
  * [58] M. Preisig, J. H. Prévost, Stabilization procedures in coupled poromechanics problems: A critical assessment, International Journal for Numerical and Analytical Methods in Geomechanics 35 (11) (2011) 1207–1225.
  * [59] L. Monforte, P. Navas, J. M. Carbonell, M. Arroyo, A. Gens, Low-order stabilized finite element for the full biot formulation in soil mechanics at finite strain, International Journal for Numerical and Analytical Methods in Geomechanics 43 (7) (2019) 1488–1515.
  * [60] C. R. Dohrmann, P. B. Bochev, A stabilized finite element method for the stokes problem based on polynomial pressure projections, International Journal for Numerical Methods in Fluids 46 (2) (2004) 183–201.
  * [61] S. Badia, E. Neiva, F. Verdugo, Linking ghost penalty and aggregated unfitted methods, Computer Methods in Applied Mechanics and Engineering 388 (2022) 114232.
  * [62] H. Kothari, R. Krause, A generalized multigrid method for solving contact problems in lagrange multiplier based unfitted finite element method, Computer Methods in Applied Mechanics and Engineering 392 (2022) 114630.
  * [63] S. Sticko, G. Ludvigsson, G. Kreiss, High-order cut finite elements for the elastic wave equation, Advances in Computational Mathematics 46 (3) (2020) 1–28.
  * [64] F. Armero, Formulation and finite element implementation of a multiplicative model of coupled poro-plasticity at finite strains under fully saturated conditions, Computer Methods in Applied Mechanics and Engineering 171 (3-4) (1999) 205–241.
  * [65] I. F. Collins, T. Hilder, A theoretical framework for constructing elastic/plastic constitutive models of triaxial tests, International journal for numerical and analytical methods in geomechanics 26 (13) (2002) 1313–1347.
  * [66] W. M. Coombs, R. S. Crouch, Algorithmic issues for three-invariant hyperplastic critical state models, Computer methods in applied mechanics and engineering 200 (25-28) (2011) 2297–2318.