#  Multiscale multimesh finite element method | $\text{M}^{2}$-FEM: Hierarchical mesh-decoupling for integral structural theories

Wei Ding∗ ∗To whom correspondence should be addressed. Email:
<EMAIL_ADDRESS>or<EMAIL_ADDRESS>Ray W. Herrick Laboratories, School
of Mechanical Engineering
Purdue University, West Lafayette, Indiana - USA, 47907  Sansit Patnaik  Ray
W. Herrick Laboratories, School of Mechanical Engineering
Purdue University, West Lafayette, Indiana - USA, 47907  Fabio Semperlotti∗
Ray W. Herrick Laboratories, School of Mechanical Engineering
Purdue University, West Lafayette, Indiana - USA, 47907

###### Abstract

This study presents a generalized multiscale multimesh finite element method
(M2-FEM) that addresses several long-standing challenges in the numerical
simulation of integral structural theories, often used to model multiscale and
nonlocal effects. The major challenges in the numerical simulation of integral
boundary value problems are primarily rooted in the coupling of the spatial
discretization of the global (parent) and integral (child) domains which
severely restricts the computational efficiency of existing algorithms by
imposing an implicit trade-off in the accuracy achieved by the child domain
and in the resources dedicated to the simulation of the overall parent domain.
One of the most defining contributions of this study consists in the
development of a mesh-decoupling technique that generates isolated sets of
meshes such that the parent and child domains can be discretized and
approximated independently. This mesh-decoupling has a multi-fold impact on
the simulation of integral theories such that, when compared to existing
state-of-the-art techniques, the proposed algorithm achieves simultaneously
better numerical accuracy and efficiency (hence allowing a greater flexibility
in both mesh size and computational cost trade-off decisions), greater ability
to adopt generalized integral kernel functions, and the ability to handle non-
regular (non-rectangular) domains via unstructured meshing. In this study, we
choose a benchmark problem based on an extended version of the Eringen’s
nonlocal elasticity theory (implicitly, a multiscale theory) that leverages
the use of generalized attenuation kernels and non-constant horizons of
nonlocality. Nonetheless, the proposed $\text{M}^{2}$-FEM algorithm is very
general and it can be applied to a variety of integral theories, even beyond
structural elasticity.

### Highlights

  * •

Finite element solver for integral structural theories: $\text{M}^{2}$-FEM

  * •

Mesh-decoupling and scale-bridging for scale-localized approximations

  * •

Ability to achieve both high accuracy and high efficiency simultaneously

  * •

Generalized kernel functions allowing modeling different material classes

  * •

Unstructured mesh to simulate irregular structures and non-similar integral
horizons

### Keywords: Integral elasticity | Multiscale effects | Nonlocal effects | Mesh-decoupling | Scale-bridging

##  Introduction

In recent years, several complex classes of heterogeneous material systems
including, for example, metamaterials [1], porous solids [2], functionally-
graded solids [3], lattice structures [4], and composites [5], have become
progressively more accessible from a fabrication perspective. These novel
classes of materials find important applications in the most diverse fields of
engineering, nanotechnology, biotechnology, and medicine. The ability to
exploit these complex material systems in practical applications requires
structural models that are capable of reliably and accurately predicting the
structural response. Several theoretical and experimental investigations have
reinforced the understanding that these materials exhibit non-negligible non-
classical behavior, that manifest at the continuum level, as multiscale
effects [5, 6, 7], nonlocal effects [8, 9, 10], and even combined multiscale-
nonlocal effects [11, 12].

The inability of the classical local continuum theory to efficiently capture
the non-classical multiscale and nonlocal effects motivated the development of
non-classical continuum theories. From a broad mathematical perspective, all
existing non-classical continuum theories can be broadly identified either as
integral or gradient continuum theories. Detailed theoretical investigations
and a diverse range of applications based over the past 50 years, have made it
increasingly evident that the integral approach is more general when compared
to the gradient approach [13, 14, 15, 16]. In the following, we provide a
high-level discussion on nonlocal and multiscale effects to briefly clarify
this latter aspect and to substantiate the choice, made later in this study,
to simulate specific integral theories.

From a physical perspective, both nonlocal and multiscale effects are
characterized by a collective and multilevel transfer of information, from the
finer structural (material and/or geometric) scales to the coarser (often,
continuum) scale. The collective information leading to nonlocal effects is
associated to multiple long-range interactions typically distributed in
parallel within a specific finer structural scale [16, 17]; in the case of
multiscale effects, they correspond to a hierarchical distribution of multiple
structural scales operating in parallel [18, 19]. The coexistence of
hierarchical structural scales and localized long-range connections at a given
scale, leads to multiscale nonlocality [11].

Integral approaches allow resolving the collective information by continuously
(in a mathematical sense) accounting the contributions via model-specific
convolution(s) over a set of interacting points localized either in a pre-
defined nonlocal horizon [15], or a set of structural scales [18], or both
simultaneously [11]. Irrespective of the specific constitutive material model,
integral approaches introduce a kernel function that characterizes the
relative decay in the multilevel exchange of information, and a horizon of
influence (as a set of material points in a localized scale, or a set of
scales, or both) that supports the kernel function. Note that, from a
mathematical standpoint, the multilevel exchange of information originating
either from multiple points (as in, nonlocal effects) and/or multiple
hierarchical scales (as in, multiscale effects), results in a similar
structure of the enriched integral constitutive equations at the continuum
scale. Consequently, in literature, nonlocal effects are often treated
implicitly as multiscale effects (and the vice versa) [18, 16, 17]. With this
understanding, in the following, we have used “multiscale effects" as a broad
term that also includes nonlocal effects.

Gradient approaches gained popularity until the past few decades mostly for
their easier numerical implementation, when compared to the significant
complexities involved in integral approaches. However, in contrast to the
seamless information resolution afforded by the integral approaches, gradient
approaches use a limited set of higher-order gradients that capture the non-
classical information from a limited number of higher-order (i.e. nearest,
next-nearest, and so on) interactions or structural scales. In fact, gradient
approaches are often derived from integral approaches by adding specialized
assumptions on the functional variation and symmetry of the kernel function
[18, 16]. These assumptions introduce major theoretical constraints in the
modeling of heterogeneous and asymmetric microstructures [13], and generalized
attenuation-dispersion behavior [14]. Consequently, gradient approaches have
found rather narrow and, on several instances, even erroneous applications
[13]. Thanks to the ability of admitting generalized kernels, integral
approaches (particularly the more recent approaches) enable the modeling of a
wide range of microstructural classes [20, 21, 22, 23] exhibiting intricate
localized (softening or stiffening) responses [11] and various attenuation-
dispersion characteristics [15, 4]. Indeed, in recent times, integral
approaches have received renewed and greater attention in non-classical
modeling.

A detailed literature survey reveals that, despite the growing attention,
numerical techniques necessary to fully realize the potential of integral
frameworks, and complement their application, are yet to be developed [7, 24].
In §1.1 below, we discuss the limitations of existing numerical techniques and
their common origin in the sharing of mesh discretization among the different
structural scales. Following this discussion, we present an overview of the
main contributions of the present study to the development of advanced
numerical techniques applicable to integral frameworks, in §1.2.

###  Limitations of existing numerical solvers for integral theories

To-date, only a handful of numerical techniques, including the NL-FEM
(nonlocal finite element method) [25, 26] and the f-FEM (fractional-order
finite element method) [27, 28], have successfully simulated generalized
integral structural theories. The most critical limitation, common to these
existing numerical strategies, consists in the explicit sharing of a globally
defined spatial mesh in the discretization of the coarse (parent) and fine
(child) length-scales (or domains). The shared discretization, or
equivalently, mesh-coupling, between different material scales adversely
impacts the simulation of multiscale problems leading to:

  * •

Inability to deliver high accuracy and efficiency, simultaneously. Mesh
coupling increasingly leads to computationally very expensive simulations with
minimal increase in accuracy. Note that the overall accuracy in multiscale
problems is strongly dependent on the spatial resolution achieved at the child
domain. Consequently, the finest spatial resolution directly determines the
size of the overall mesh, as a result of the underlying spatial coupling.
Hence, in applications involving a greater disparity between coarse and fine
scales, this restriction could lead to a very computationally expensive mesh
at the level of the coarse domain, potentially even preventing the ability to
numerically solve the problem. As an example, consider a porous beam with very
fine pores’ scales compared to the characteristic length of the beam. In this
case, it is not difficult to envision how a fine discretization for the child
domain could lead to a multi-fold increase in the mesh size (in terms of the
number of elements) at the level of the parent (i.e. beam scale) domain,
drastically impacting the overall efficiency of the simulation technique.

  * •

Inability to simulate irregular structural configurations. The direct sharing
of the spatial discretization implicitly requires the parent and child domains
to be geometrically similar which, in practice, restricts the application of
existing approaches to the simulation of multiscale structures with either
rectangular domains (such as, for example, beams and plates) or at best,
structures that can be reduced parametrically to rectangular domains (such as
rectangular panels and shells). Consequently, it is not surprising that the
application of existing approaches to multiscale structures with irregular
material geometry and (geometrically) non-similar integral (wither multiscale
or nonlocal) horizon geometries are absent from the literature.

Notably, in addition to the mesh-coupling, existing approaches are not well
equipped to directly accommodate the most generalized (singular) kernel
functions such as, for example, the power-law kernel (in constant- [15],
variable- [29], or distributed-order [11] version) which has found widespread
application in the recent times [20, 21, 22]. f-FEM, that simulates differ-
integral partial differential equations based on fractional
calculus111Fractional calculus is a branch of mathematics that enables the
differentiation and integration of functions beyond the classical set of
natural numbers and extends it to real and complex numbers, and even
functions. A discussion of fractional-order derivatives and integrals is
beyond the scope of this study and the interested reader is referred to [30]
for detailed discussions., does admit singular kernels, but it still adopts
coupled meshing.

The several limitations discussed above largely restrict the application of
existing computational strategies to the simulation of highly multiscale
structures characterized by strong underlying heterogeneity [31, 18]. An
analysis of the most recent developments in the foundational theory and
application of multiscale models reveals that classical phenomenological
approaches (that can be more readily simulated by existing solvers) are of
limited practical significance in the several novel classes of materials
mentioned previously. The need for accuracy has suggested embracing the
multiscale heterogeneous nature of these materials [20, 12, 32], as opposed to
a more classical homogeneous and symmetric view typically associated with
phenomenological (representative volume element, RVE-based) approaches. Some
notable studies that provide concrete strategies for this renewed effort
include mesoscale approaches informed by extensive molecular dynamics
simulations [19, 18], nonlocal approaches based on variable-order [29] and
distributed-order [11] fractional calculus, and even, data-driven methods [33,
34].

The large material parameter space associated with the aforementioned advanced
theories ultimately necessitates the application of advanced inverse
techniques (e.g. machine learning-based techniques) for constitutive discovery
[35, 33, 29, 34]. The extensive datasets used in these discovery processes,
are often generated artificially by the associated numerical solvers based
either on a priori knowledge of selective parameter combinations [33, 36] or
via novel strategies to implement one-shot discovery [20, 34]. To this end,
the prohibitive computational cost required for better accuracy, the limited
ability in modeling various geometries resulting from the mesh-coupling, and
the inability in admitting generalized kernel functions, drastically impact
the generation of critical training datasets and in general, even the forward
implementation of these advanced theories in future applications. Evidently,
existing multiscale computational approaches are not well-suited for the next-
generation of applications.

###  Overall approach and contributions of the present study

The previous discussion motivates the need to revisit the mesh topology and
sharing strategy adopted in existing solvers. It is evident that the ideal
computational approach, which addresses the limitations outlined in §1.1,
requires an effective decoupling between the spatial discretization of the
parent and child domains. It is important to understand that while the
decoupling takes place in the numerical discretization, from a theoretical
standpoint the information exchange across scales is still maintained. This
condition can be achieved independently of the meshing strategy by ensuring
that the adopted multiscale theory is thermodynamically consistent [37, 18],
and by imposing appropriate connectivity among co-existing multiscale sets of
material points in the numerical simulation (more commonly known as scale-
bridging) [24, 38]. These concepts define the broad approach of the present
study.

One of the most significant contribution of this study, towards advanced
integral simulation, is the development of a novel mesh-decoupling technique
that lays the foundation of a multiscale multimesh finite element method,
denoted as $\text{M}^{2}\text{-FEM}$.

Recall that the weak-form approximation of integral theories consists of the
approximation of (at least) two nested integrals, the parent-domain integral
and the child-domain constitutive (multiscale/nonlocal) integral. The proposed
mesh-decoupling technique generates isolated sets of meshes in the parent and
child domains. By approximating the nested constitutive integral in the
decoupled child mesh, the integral constitutive information is encapsulated
within the child domain. This approach minimizes the transmission of numerical
error generated from the child domain approximation to the parent domain (and
likewise in the reverse direction). Thus, mesh-decoupling ultimately decouples
the errors in the approximation of the two domain-specific integrals. The most
significant outcome of this decoupling is the remarkable ability of
$\text{M}^{2}$-FEM to independently tune parent and child mesh sizes, to
achieve the desired efficiency and accuracy, simultaneously. Further,
decoupling also allows unstructured meshing of both domains to model irregular
configurations, and tailored numerical algorithms to embrace generalized
singular kernels while approximating child domain integrals. In the present
work, the development of $\text{M}^{2}\text{-FEM}$ is complemented by the
following additional set of studies:

  * •

Validation and convergence analyses are performed through a series of targeted
case studies to validate the implementation of the mesh-decoupling mechanism
and establish convergence conditions. The role of the child and parent domain
resolutions, and the kernel function in determining the accuracy and
convergence are investigated. As we will show, $\text{M}^{2}$-FEM enables a
rapidly convergent and accurate simulation with relatively coarse parent
domain mesh; a remarkable feature that is not attainable with existing
solvers.

  * •

Computational complexity analysis and simulation of non-rectangular structures
are presented to showcase the potential impact of $\text{M}^{2}$-FEM on the
next-generation of multiscale problems. We will perform specific case studies
to demonstrate how the superior convergence characteristics enable a relative
tuning of the parent and child mesh resolutions to achieve superior accuracy
and faster simulation times, simultaneously. Finally, a series of advanced
geometries, identified from practically relevant structures in a variety of
fields, are considered to showcase the abilities of $\text{M}^{2}$-FEM;
existing solvers cannot simulate the considered geometries.

For all the simulations presented in this study, we have adopted as the
underlying structural theory an extended version of Eringen’s integral
(strain-driven) nonlocal theory [17] that admits non-singular and singular
kernel functions. Note that Eringen’s nonlocal theory is implicitly multiscale
in nature [18]. To this end, we emphasize that the algorithm developed is very
general and can be directly applied to any integral multiscale formulation
related to structural mechanics or other continuum theories. We decided to
adopt Eringen’s approach because the broader scientific community still
relates more closely to this seminal theory. It is also important to note that
critical aspects, such as model convexity, thermodynamic consistency, and
mathematical well-posedness are difficult to address when using the strain-
driven approach [16]. We merely note that the recently proposed displacement-
driven approach [14] is a more suitable candidate since it is intrinsically
convex, thermodynamically consistent, and well-posed. Nonetheless, since model
convexity is essential for finite element solvers, the successful simulation
of the strain-driven theory in this study, essentially guarantees a rather
straightforward application of $\text{M}^{2}$-FEM also to the simulation of
the displacement-driven theory.

The remainder of this paper is structured as follows. The mesh-decoupling
technique and the scale-bridging algorithms are developed in §2, alongside the
finite element approximations of both parent and child domains. Next, in §3,
we develop the matrix representation of the approximated components, assemble
the domain-specific approximations, and derive the algebraic expressions for
computer implementation. Results of validation and convergence studies are
presented in §4. Simulations obtained from a series of practically relevant
multiscale configurations are presented in §5 to highlight the envisioned
impact of $\text{M}^{2}$-FEM. Finally, a brief summary and the future outlook
are presented in §6.

##  Mathematical formulation: decoupled multiscale discretization and
numerical approximation

In this section, we develop the mesh-decoupling algorithm that is at the
foundation of M2-FEM. This development includes the definition of separate
meshes for the discretization of different material scales as well as the use
of the different meshes to derive the multiscale finite element approximations
for the simulation of the background theory. Note that, while the algorithm is
general, we will derive the approximations for a 2D problem (as opposed to the
full 3D problem) to facilitate the understanding of important aspects of the
algorithm. Here below, we briefly summarize Eringen’s integral theory
(focusing on the constitutive relations and governing equations) and the
notation scheme adopted in the study, to more efficiently utilize them in
deriving and discussing the discretization and approximation techniques.

###  Brief summary of the adopted integral theory

According to the strain-driven approach to nonlocal elasticity[17], the
stress-strain constitutive relation at any material point in a nonlocal solid,
can be expressed using integral operators in the following form:

$\bm{\sigma}=\bm{C}:\underbrace{\int_{\mathcal{H}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\bm{\epsilon}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}}_{\textrm{nonlocal
integral}}$ (1)

where $\bm{\sigma}$, $\bm{C}$, and $\bm{\epsilon}$ are the nonlocal stress,
stiffness, and local strain tensors, respectively. Further, $\bm{x}$ denotes a
specific material point on the nonlocal solid (say, $\Omega$); and
$\bm{x}^{\prime}\in\mathcal{H}$ is a dummy spatial variable denoting the
(remaining) materials points that interact with $\bm{x}$ via long-range
connections and are contained within a characteristic nonlocal horizon
$\mathcal{H}$. The strength of the nonlocal effects are captured by the
nonlocal kernel function denoted as $\mathcal{K}(\bm{x},\bm{x}^{\prime})$ with
support on the nonlocal horizon $\mathcal{H}$.

The strong form of the governing equations can be derived by variational
minimization of the potential energy functional obtained from the constitutive
framework presented in Eq. (1), as [16]:

$\displaystyle\bm{\nabla}\cdot\bm{\sigma}+\rho\bm{f}-\rho\bm{\ddot{u}}$
$\displaystyle=0,\quad\forall\bm{x}\in\bm{\Omega}$ (2)

The governing equations are subject to the following boundary and initial
conditions:

$\displaystyle\delta\bm{u}=0\quad\text{or}\quad\bm{\sigma}\cdot\bm{n}-\bm{t}=0,\quad\forall\bm{x}\in\bm{\Gamma}$
(3)
$\displaystyle\delta\bm{u}=0\quad\text{and}\quad\delta\dot{\bm{u}}=0,\quad\forall\bm{x}(t=0)$

where $\bm{\Gamma}$ denotes the boundary of the nonlocal solid, and $\bm{t}$
denotes the traction force acting on the surface whose outward normal is
indicated by $\bm{n}$. Analogous to the classical approach, we use the weak
form representation of Eq. (2), presented below (see [16] for derivation), for
finite element discretization and approximation:

$\displaystyle\int_{\bm{\Omega}}\bm{\sigma}:\bm{\nabla}\delta\bm{u}\mathrm{d}\bm{\Omega}-\int_{\bm{\Omega}}\rho(\bm{f}-\ddot{\bm{u}})\delta\bm{u}\mathrm{d}\bm{\Omega}-\int_{\partial\bm{\Omega}}\bm{\sigma}\cdot\bm{n}\delta\bm{u}\mathrm{d}\bm{\Gamma}$
$\displaystyle=0$ (4)

In proceeding further, we assume that the nonlocal kernel function
$\mathcal{K}(\bm{x},\bm{x}^{\prime})$ is positive-definite and symmetric such
that the total strain energy functional is convex and the governing equations
in Eq. (2) are well-posed [16].

It is evident from the Eqs. (1,4) that the finite element implementation of
the assumed background theory requires the approximations of different
physical variables defined on separate material scales. These approximations
correspond either to: 1) the weak-form integral over the material domain
($\bm{\Omega}$), and 2) the convolution integral defined over the nonlocal
horizon ($\mathcal{H}$). Henceforth, we simply refer to the two material
scales as the parent domain (that is, the material domain) and the child
domain (that is, the nonlocal horizon).

###  Remarks on the adopted notation scheme

The separate meshing schemes used to discretize the two different domains and
the integral nature of the underlying formulation result in a large number of
variables and finite element parameters. In the interest of streamlining the
following derivations and discussions, as well as to facilitate the reading,
we discuss the notation scheme here below. The most general notation for a
given variable, say $\square$, that encompasses the subsequent development can
be expressed as:

$\prescript{\mathrm{u_{-}}}{l_{-}}{}\square^{\mathrm{u_{+}}}_{l_{+}}$ (5)

where the sets of superscripts $\\{\mathrm{u}_{-},\mathrm{u}_{+}\\}$ and
subscripts $\\{l_{-},l_{+}\\}$ can be specified to denote different variables.

We typeset the superscripts in Roman font and use them to generate symbols for
different variables well-defined either on the decoupled-meshes or on the
solid domain. The complete list of superscripts is here below:

  * •

$\\{\mathrm{p},\mathrm{c}\\}$ represent quantities defined on the parent and
child domain.

  * •

$\\{\mathrm{b},\mathrm{u},\mathrm{t}\\}$ represent quantities defined on the
complete boundary, and the segments of the boundary subjected to displacement
and traction loads.

For enhanced readability, we typeset the subscript using italic fonts and use
them for indexing. Besides the commonly used indices $i$ and $j$, we define
the following indices to represent specific physical quantities:

  * •

$\\{r,\tilde{r},s\\}$ are used for parent, dummy parent (defined in
$\mathcal{H}$, see §2.4.2), and child element indexing.

  * •

$\\{d,g\\}$ are used for nodal and Gauss point indexing.

  * •

$\\{u,v\\}$ are used for number of degrees of freedom (DOF) and weighting
function indexing.

Further, we use the right-hand scripts $\\{\mathrm{u}_{+},l_{+}\\}$ to
indicate quantities defined at the current scale (for example,
$\Omega_{r}^{\textrm{p}}$ and $\Omega_{s}^{\textrm{c}}$ that define finite
elements in the parent domain and child domain, respectively) and left-hand
scripts $\\{\mathrm{u}_{-},l_{-}\\}$ to indicate coarser scale quantities.
This convention proves effective when describing multiscale quantities.
Consider for example, $\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$,
which represents finite elements in the child domain (see §2.4.1 or Eq. (22)).
In $\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$, the right scripts
${\square}_{s}^{\textrm{c}}$ specify that
$\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$ is the $s^{\mathrm{th}}$
element in the child domain, and the left scripts
$\prescript{\textrm{p}}{rg}{\square}$ specify that this child element is
defined at the nonlocal horizon of $g^{\mathrm{th}}$ parent Gauss point in the
$r^{\mathrm{th}}$ parent element.

Two final remarks on the notation: 1) we have used a comma to denote partial
derivatives in agreement with the classical tensor algebra notation [39]; and
2) in the interest of a general derivation, we have used different symbols,
$\Phi(\cdot)$ and $\Psi(\cdot)$, to denote the shape functions used to
interpolate the geometrical coordinates and the displacement degrees of
freedom, respectively; the remaining variables are introduced as they appear
in the derivation below.

###  Parent domain discretization and approximation

We first develop the finite element approximation of the weak-form governing
equations over the parent domain $\bm{\Omega}$ and boundary $\bm{\Gamma}$.
Note that, since the governing equations on the material scale presented in
Eq. (2) exhibit a functional character similar to their classical (local)
counterpart, the parent domain approximations very closely follows a classical
finite element approximation of the 2D plane strain problem. However, we fully
derive the parent domain approximations, following a discretization of the
weak-form integral at the parent level in order to: 1) clearly highlight the
mesh-decoupling in action and the connectivity algorithm developed to preserve
the multiscale exchange of information, and 2) enable an efficient discussion
on the matrix representation and integration of the dual-level approximations
in the final implementation in §3.

####  Discretization

Consider the nonlocal elastic material occupying, in its undeformed state, a
2D Euclidean parent domain $\bm{\Omega}^{\textrm{p}}$ referred to orthogonal
Cartesian coordinates $\bm{x}=(x,y)$. The 2D domain is discretized into
$N^{\textrm{p}}$ disjoint finite elements $\Omega_{r}^{\textrm{p}}$ such that:

$\bigcup_{r=1}^{N^{\textrm{p}}}\Omega_{r}=\bm{\Omega},\quad\Omega_{i}\cap\Omega_{j}=\varnothing,~{}\forall
i\neq j$ (6)

where the superscript p indicates quantities related to the parent domain. In
all the computational developments considered in the following, the parent
mesh is assumed to be non-rectangular in shape, unless specified otherwise.
The boundary $\bm{\Gamma}$ is also discretized into $N^{\textrm{b}}$ elements
that can be differentiated into the following different types depending on the
nature of the boundary condition:

$\displaystyle\bigcup_{r=1}^{N^{\textrm{b}}}\Gamma_{r}$
$\displaystyle=\bigcup_{i=1}^{N^{\textrm{u}}}\Gamma_{i}^{\textrm{u}}+\bigcup_{j=1}^{N^{\textrm{t}}}\Gamma_{j}^{\textrm{t}}=\bm{\Gamma},\quad$
$\displaystyle\Gamma_{i}^{\textrm{u}}\cap\Gamma_{j}^{\textrm{t}}=\varnothing,~{}\forall
i\in[1,N^{\textrm{u}}],~{}j\in[1,N^{\textrm{t}}]$ (7)
$\displaystyle\bigcup_{r=1}^{N^{\textrm{u}}}\Gamma_{r}^{\textrm{u}}$
$\displaystyle=\bm{\Gamma}^{\textrm{u}},\quad$
$\displaystyle\Gamma_{i}^{\textrm{u}}\cap\Gamma_{j}^{\textrm{u}}=\varnothing,~{}\forall
i\neq j$ $\displaystyle\bigcup_{r=1}^{N^{\textrm{t}}}\Gamma_{r}^{\textrm{t}}$
$\displaystyle=\bm{\Gamma}^{\textrm{t}},\quad$
$\displaystyle\Gamma_{i}^{\textrm{t}}\cap\Gamma_{j}^{\textrm{t}}=\varnothing,~{}\forall
i\neq j$

where $N^{\textrm{u}}$ and $N^{\textrm{t}}$ are the number of boundaries
subjected to the displacement boundary condition (represented by
$\bm{\Gamma}^{\textrm{u}}$) and traction boundary conditions (represented by
$\bm{\Gamma}^{\textrm{t}}$), respectively.

Using the above parent domain discretization, the continuous integral
corresponding to the simplified variational governing equations (see Eq. (4))
is approximated as the limit of the following summation:

$\sum_{r=1}^{N^{\textrm{p}}}\int_{\Omega_{r}^{\textrm{p}}}\bm{\sigma}:\bm{\nabla}_{\bm{x}}\delta\bm{v}\mathrm{d}\bm{\Omega}-\sum_{r=1}^{N^{\textrm{p}}}\int_{\Omega_{r}^{\textrm{p}}}\rho(\bm{f}-\ddot{\bm{u}})\cdot\delta\bm{v}\mathrm{d}\bm{\Omega}-\sum_{r=1}^{N^{\textrm{t}}}\int_{\Gamma_{r}^{\textrm{t}}}\bm{\sigma}\cdot\bm{n}\delta\bm{v}\mathrm{d}\bm{\Gamma}=0$
(8)

where $\bm{u}$ and $\delta\bm{v}$ denote the nodal degrees of freedom (DOF)
and the weighting function. Further, the boundary conditions are expressed as:

$\displaystyle\bm{u}|_{\bm{\Gamma}^{\textrm{u}}}$
$\displaystyle=\overline{\bm{u}}$ (9)
$\displaystyle\bm{\sigma}\cdot\bm{n}|_{\bm{\Gamma}^{\textrm{t}}}$
$\displaystyle=\overline{\bm{t}}$

where $\overline{\bm{u}}$ and $\overline{\bm{t}}$ denote the externally
applied displacement and traction loads. By combining the approximations
presented above for the governing equations and boundary conditions, the first
step discretization developed here at the parent domain
$\bm{\Omega}^{\textrm{p}}$ results in the following approximation of the
governing variational equation:

$\displaystyle\sum_{r=1}^{N^{\textrm{p}}}\int_{\Omega_{r}^{\textrm{p}}}\bm{\sigma}:\bm{\nabla}_{\bm{x}}\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}-\sum_{r=1}^{N^{\textrm{p}}}\int_{\Omega_{r}^{\textrm{p}}}\rho(\bm{f}-\ddot{\bm{u}})\cdot\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}-\sum_{r=1}^{N^{\textrm{t}}}\int_{\Gamma_{r}^{\textrm{t}}}\overline{\bm{t}}\cdot\delta\bm{v}\mathrm{d}\Gamma_{r}^{\textrm{t}}$
$\displaystyle=0$ (10)

####  Numerical integration

As evident from Eq. (10), the approximation of the variational equation
consists of two different components: 1) the numerical evaluation of the
discretized integral expression over the parent domain elements
$\Omega_{r}^{\textrm{p}}$, and 2) the numerical evaluation of the discretized
integral expression over the boundary elements $\Gamma_{r}^{\textrm{t}}$. In
the following, we first evaluate the parent domain integrals and then consider
the boundary domain integrals.

In order to perform the numerical integration, we adopt an isoparametric
formulation and introduce a 2D natural coordinate system. In more specific
terms, for a given parent element $\Omega_{r}^{\textrm{p}}$, $\forall
r\in[1,N^{\textrm{p}}]$, we transform the material points $\bm{x}=(x,y)$ in
$\Omega_{r}^{\textrm{p}}$ from the 2D Cartesian coordinates to natural
coordinates
$\bm{\xi}_{r}^{\textrm{p}}=(\xi_{r}^{\textrm{p}},\eta_{r}^{\textrm{p}})$. This
isoparametric mapping can be expressed as:

$\bm{x}(\bm{\xi}_{r}^{\textrm{p}})=\sum_{d=1}^{d^{\textrm{p}}}\Phi_{d}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\bm{x}_{rd}^{\textrm{p}}=\\{\Phi_{d}(\bm{\xi}_{r}^{\textrm{p}})\\}^{\textrm{T}}\\{\bm{x}_{rd}^{\textrm{p}}\\}$
(11)

such that the parent element $\Omega_{r}^{\textrm{p}}$ can be transformed into
a rectangular element $\hat{\Omega}_{r}^{\textrm{p}}$ in the natural
coordinate system. In the latter transformation,
$\Phi_{d}^{\textrm{p}}(\cdot)$ and $d^{\textrm{p}}$ are the shape functions
and the number of nodes in the parent element $\Omega_{r}^{\textrm{p}}$,
respectively. Finally, $\\{\Phi_{d}(\bm{\xi}_{r}^{\textrm{p}})\\}$ and
$\\{\bm{x}_{rd}^{\textrm{p}}\\}$ (where $d=1,2,...,d^{\textrm{p}}$), that
denote the vectors assembling the nodal shape functions and the corresponding
nodal coordinates, respectively, are defined as:

$\displaystyle\\{\Phi_{d}(\bm{\xi}_{r}^{\textrm{p}})\\}$
$\displaystyle=\\{\Phi_{1}(\bm{\xi}_{r}^{\textrm{p}}),\Phi_{2}(\bm{\xi}_{r}^{\textrm{p}}),...,\Phi_{d^{\textrm{p}}}(\bm{\xi}_{r}^{\textrm{p}})\\}^{\textrm{T}}$
(12a) $\displaystyle\\{\bm{x}_{rd}^{\textrm{p}}\\}$
$\displaystyle=\left\\{\\{x_{r}^{\textrm{p}},y_{r}^{\textrm{p}}\\}_{1}^{\textrm{T}},\\{x_{r}^{\textrm{p}},y_{r}^{\textrm{p}}\\}_{2}^{\textrm{T}},...,\\{x_{r}^{\textrm{p}},y_{r}^{\textrm{p}}\\}_{d^{\textrm{p}}}^{\textrm{T}}\right\\}^{\textrm{T}}$
(12b)

where the subscript $r$ denote quantities supported on the $r^{\textrm{th}}$
parent element. Note that, analogous to classical finite element methodology,
the isoparametric transformation between the Cartesian and the natural
coordinate systems is affine in nature [39]. More specifically, the shape
functions $\Phi_{d}^{\textrm{p}}(\cdot)$ used for geometrical interpolation
should be linear functions when defined in either coordinate systems; see
additional details in §2.4.2. We consider Gauss points
$\bm{\xi}_{rg}^{\textrm{p}}=(\xi_{rg}^{\textrm{p}},\eta_{rg}^{\textrm{p}})$
and the corresponding weights $w_{rg}^{\textrm{p}}$ in the isoparametric
element $\hat{\Omega}_{r}^{\textrm{p}}$. Now by using the classical Gauss-
Legendre quadrature principles, the discretized integrals in Eq. (10) over the
parent element $\Omega_{r}^{\textrm{p}}$, with Gauss points
$\bm{\xi}_{rg}^{\textrm{p}}=(\xi_{rg}^{\textrm{p}},\eta_{rg}^{\textrm{p}})$
and the corresponding weights $w_{rg}^{\textrm{p}}$, can be approximated as:

$\displaystyle\int_{\Omega_{r}^{\textrm{p}}}\bm{\sigma}:\bm{\nabla}_{\bm{x}}\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}$
$\displaystyle=\sum_{g=1}^{g^{\textrm{p}}}{w}_{rg}^{\textrm{p}}|\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|\bm{\sigma}(\bm{\xi}_{rg}^{\textrm{p}}):\bm{\nabla}_{\bm{x}}\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})$
(13a)
$\displaystyle\int_{\Omega_{r}^{\textrm{p}}}\rho(\bm{f}-\ddot{\bm{u}})\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}$
$\displaystyle=\sum_{g=1}^{g^{\textrm{p}}}w_{rg}^{\textrm{p}}|\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|\rho(\bm{\xi}_{rg}^{\textrm{p}})(\bm{f}(\bm{\xi}_{rg}^{\textrm{p}})-\ddot{\bm{u}}(\bm{\xi}_{rg}^{\textrm{p}}))\cdot\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})$
(13b)

where $g^{\textrm{p}}$ denotes the number of Gauss points.
$|\bm{J}_{g}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|$ denotes the Jacobian
of the transformation at the Gauss point $\bm{\xi}_{rg}^{\textrm{p}}$ and is
evaluated in the following manner:

$\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})=\frac{\partial\bm{x}(\bm{\xi}_{rg}^{\textrm{p}})}{\partial\bm{\xi}_{r}^{\textrm{p}}}=\sum_{d=1}^{d^{\textrm{p}}}\frac{\partial\Phi_{d}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})}{\partial\bm{\xi}_{r}^{\textrm{p}}}\bm{x}_{rd}^{\textrm{p}}=\sum_{d=1}^{d^{\textrm{p}}}\bm{x}_{rd}^{\textrm{p}}\otimes\Phi_{d,\bm{\xi}_{r}^{\textrm{p}}}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})=\\{\Phi_{d,\bm{\xi}_{r}^{\textrm{p}}}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{\textrm{T}}\\{\bm{x}_{rd}^{\textrm{p}}\\}$
(14)

in which $\otimes$ denotes the outer product between two tensors (note that
the subscripts are indices of $\bm{x}$ and $\Phi$, not indices of components
of these tensors). As we have clarified in §2.2, the subscript
$\bm{\xi}_{r}^{\textrm{p}}$ (that consists of $\xi_{r}^{\textrm{p}}$ and
$\eta_{r}^{\textrm{p}}$) indicates the partial derivative with respective to
the isoparametric variable in the $r^{\textrm{th}}$ parent domain.

Note that the expressions derived in Eq. (13) require an evaluation of the
displacement ($\ddot{\bm{u}}$) and the weighting function ($\delta\bm{v}$ and
$\bm{\nabla}_{\bm{x}}\delta\bm{v}$). By adopting the classical Galerkin
approach, the weighting functions are chosen as the shape functions adopted to
interpolate the displacement variables. Consequently, the displacement and
weighting functions can be interpolated in terms of their corresponding nodal
quantities in the following manner:

$\displaystyle\ddot{\bm{u}}(\bm{\xi}_{rg}^{\textrm{p}})$
$\displaystyle=\sum_{u=1}^{m^{\textrm{p}}}\Psi_{u}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\ddot{\bm{u}}_{ru}^{\textrm{p}}$
(15a) $\displaystyle\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})$
$\displaystyle=\sum_{v=1}^{m^{\textrm{p}}}\Psi_{v}^{p}(\bm{\xi}_{r}^{p})\delta\bm{v}_{rv}^{\textrm{p}}$
(15b)

where subscripts $u,v\in\\{1,2,...,m^{\textrm{p}}\\}$ are the two indices for
$\bm{u}$ and $\delta\bm{v}$, respectively; and $m^{\textrm{p}}$ is the total
number of DOF within the parent element. Now, by using the shape function
$\Psi(\cdot)$, the gradient of the weighting functions
$\bm{\nabla}_{\bm{x}}\delta\bm{v}$ can be immediately approximated as:

$\bm{\nabla}_{\bm{x}}\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})=\sum_{v=1}^{m^{\textrm{p}}}\delta\bm{v}_{rv}\otimes\bm{\nabla}_{\bm{\xi}_{r}^{\textrm{p}}}\Psi_{v}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\frac{\partial\bm{\xi}_{rg}^{\textrm{p}}}{\partial\bm{x}}=\sum_{v=1}^{d^{\textrm{p}}}\delta\bm{v}_{rv}\otimes\left[\bm{\nabla}_{\bm{\xi}_{r}^{\textrm{p}}}\Psi_{v}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})(\bm{J}_{rg}^{\textrm{p}})^{-1}\right]$
(16)

where the italic letters $u$ and $v$ in the subscripts serve as indices to
distinguish between shape functions for solutions (such as the term
$\ddot{\bm{u}}_{ru}^{\textrm{p}}$) and weight function terms (such as the term
$\delta\bm{v}_{rv}^{\textrm{p}}$), respectively. Substituting Eqs. (15,16)
within Eq. (13), we finally obtain the approximation of the integral over
$\Omega_{r}^{\textrm{p}}$ as:

$\displaystyle\int_{\Omega_{r}^{\textrm{p}}}\bm{\sigma}:\bm{\nabla}_{\bm{x}}\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}$
$\displaystyle=\sum_{g=1}^{g^{\textrm{p}}}{w}_{rg}^{\textrm{p}}|\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|\bm{\sigma}(\bm{\xi}_{rg}^{\textrm{p}}):\left[\sum_{v=1}^{d^{\textrm{p}}}\delta\bm{v}_{rv}\otimes\left[\bm{\nabla}_{\bm{\xi}_{r}^{\textrm{p}}}\Psi_{v}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})(\bm{J}_{rg}^{\textrm{p}})^{-1}\right]\right]$
(17a)
$\displaystyle\int_{\Omega_{r}^{\textrm{p}}}\rho(\bm{f}-\ddot{\bm{u}})\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}$
$\displaystyle=\sum_{g=1}^{g^{\textrm{p}}}w_{rg}^{\textrm{p}}|\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|\rho(\bm{\xi}_{rg}^{\textrm{p}})\left[\bm{f}(\bm{\xi}_{rg}^{\textrm{p}})-\sum_{u=1}^{d^{\textrm{p}}}\Psi_{u}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\ddot{\bm{u}}_{ru}^{\textrm{p}}\right]\cdot\left[\sum_{v=1}^{d^{\textrm{p}}}\Psi_{v}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\delta\bm{v}_{rv}^{\textrm{p}}\right]$
(17b)

By retracing the procedure outlined above, the integrals over the boundary
elements in Eq. (10) are evaluated numerically, here below. For this purpose,
the line elements in the Cartesian system (of the boundary) are mapped on to a
1D natural coordinate system via the following affine transformation:

$\bm{x}(\bm{\xi}_{r}^{\textrm{t}})=\sum_{d=1}^{d^{\textrm{t}}}\Phi_{d}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}})\bm{x}_{rd}^{\textrm{t}}$
(18)

where $\Phi_{d}^{\textrm{t}}(\cdot)$ and $d^{\textrm{t}}$ are the $C^{0}$
Lagrangian shape function and the number of nodes in the boundary element
$\Omega_{r}^{\textrm{t}}$, respectively. Again, by using Gauss-Legendre
quadrature in isoparametric coordinates, the integral over the parent element
$\Gamma_{r}^{\textrm{t}}$ shown in Eq. (10) can be approximated as:

$\int_{\Gamma_{r}^{\textrm{t}}}\overline{\bm{t}}\delta\bm{v}\mathrm{d}\Gamma_{r}^{t}=\sum_{g=1}^{g^{\textrm{t}}}w_{rg}^{\textrm{t}}|\bm{J}_{r}^{\textrm{t}}|\overline{\bm{t}}(\bm{\xi}_{rg}^{\textrm{t}})\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{t}})$
(19)

where $g^{\textrm{t}}$ denotes the number of Gauss points at
$\bm{\xi}_{rg}^{\textrm{t}}$ and the corresponding Gaussian weights
$w_{rg}^{\textrm{t}}$. The Jacobian of the transformation, denoted above as
$|\bm{J}_{r}^{\textrm{t}}|$, is obtained as:

$\bm{J}_{r}^{\textrm{t}}=\frac{\partial\bm{x}}{\partial\bm{\xi}_{r}^{\textrm{t}}}=\sum_{d=1}^{d^{\textrm{t}}}\frac{\partial\Phi_{d}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}})}{\partial\bm{\xi}_{r}^{\textrm{t}}}\bm{x}_{rd}^{\textrm{t}}=\sum_{d=1}^{d^{\textrm{t}}}\bm{x}_{rd}^{\textrm{t}}\otimes\Phi_{d,\bm{\xi}_{r}^{\textrm{t}}}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}})$
(20)

Finally, by using the approximation of the weighting function term
$\delta\bm{v}$ from Eq. (15), the boundary integrals over
$\Gamma_{r}^{\textrm{t}}$ are evaluated as:

$\int_{\Gamma_{r}^{\textrm{t}}}\overline{\bm{t}}\delta\bm{v}\mathrm{d}\Gamma_{r}^{\textrm{t}}=\sum_{g=1}^{g^{\textrm{t}}}w_{rg}^{\textrm{t}}|\bm{J}_{r}^{\textrm{t}}|\overline{\bm{t}}(\bm{\xi}_{rg}^{\textrm{t}})\left[\sum_{v=1}^{m^{\textrm{t}}}\Psi_{v}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}})\delta\bm{v}_{rv}^{\textrm{t}}\right]$
(21)

This concludes the finite element approximation at the parent domain level.

##### Resolution from parent to child domain.

Note that the evaluation of the stress $\bm{\sigma}$ requires the finite
element approximation of the integral constitutive relation in Eq. (1). The
latter convolution integral requires an additional level of approximation that
corresponds to the nonlocal horizon, that is, the child domain. This
understanding brings us to the next step, that is the finite element
discretization and approximations within the child domain.

###  Child domain discretization and approximation

In this section, we focus on the evaluation of nonlocal stress $\bm{\sigma}$
at each Gauss point in parent elements $\Omega_{rg}^{\textrm{p}}$ and boundary
elements $\Gamma_{r}^{\textrm{p}}$. As evident from the integral constitutive
relation in Eq. (1), the numerical evaluation of the stress requires the
information on the strain level at material points located within the nonlocal
horizon $\mathcal{H}$. The evaluation of this nested integral poses several
challenges when adopting an independent (child) mesh, decoupled from the
parent mesh. Some of these challenges include the approximation of generalized
(particularly, non-singular) nonlocal kernel functions and enforcing well-
defined multiscale connectivity with the parent domain.

Before developing the discretization for the child mesh, we emphasize that the
discretization and approximation of the child domain are directly influenced
by the degree of singularity of the nonlocal kernel. Non-singular kernels (for
example, exponential kernels, typical of classical nonlocal approaches) lead
to proper convolution integrals over the nonlocal horizon $\mathcal{H}$ that
can be accurately approximated using the Gauss-Legendre quadrature principle.
On the other hand, weakly-singular kernels (for example, power law kernels
with exponents $\alpha\in(0,1)$, typical of fractional calculus based methods)
lead to improper integrals that require the Gauss-Jacobi quadrature principle
along with restrictions on the location of the quadrature points [40].
Similarly, more advanced quadrature schemes could be used for strongly-
singular nonlocal kernels (for example, power law kernels with exponents
$\alpha>1$), these kernels are yet to find meaningful applications except in
applications concerned with plastic damage such as, cracks and dislocations
[18, 19]. In this study, we only consider elastic deformation across the
scales and restrict the analysis to non-singular and weakly-singular nonlocal
kernels.

####  Discretization

We consider the $r^{\mathrm{th}}$ parent element $\Omega_{r}^{\textrm{p}}$,
the discretization of nonlocal horizon for each Gauss point
$\bm{x}(\bm{\xi}_{rg}^{\textrm{p}})$, $g={1,2,...,g^{\mathrm{p}}}$ in
$\Omega_{r}^{\textrm{p}}$, consists of the following two steps in the outlined
sequence:

  * •

Step 1: determine the nonlocal horizon, that is the child domain
$\prescript{\mathrm{p}}{rg}{}\bm{\Omega}^{\mathrm{c}}=\mathcal{H}(\bm{x}(\bm{\xi}_{rg}^{\mathrm{p}}))$,
at the Gauss point. The prescript $\prescript{\textrm{p}}{rg}{}(\cdot)$
indicates quantities that belong to the $g^{\mathrm{th}}$ Gauss point in the
$r^{\mathrm{th}}$ parent element.

  * •

Step 2: discretize the nonlocal horizon to obtain the child mesh.

Step 1 above determines the shape and geometry of the child domain
$\prescript{\textrm{p}}{rg}{}\bm{\Omega}^{\textrm{c}}$ in the Euclidean space.
Since the nonlocal horizon is contained within the material boundary
$\bm{\Gamma}$, $\prescript{\textrm{p}}{rg}{}\bm{\Omega}^{\textrm{c}}$ at
different Gauss points (particularly for Gauss points
$\bm{x}(\bm{\xi}_{rg}^{\textrm{p}})$ near the material boundary) are
truncated. Figure (4-b), which shows the mesh configuration of examples (which
are simulated in §4.1), illustrates how the material boundary affects the
shape of nonlocal horizon.

Similar to parent domain, the 2D child domain
$\prescript{\textrm{p}}{rg}{}\bm{\Omega}^{\textrm{c}}$ is discretized into
$\prescript{\textrm{p}}{rg}{}N^{\textrm{c}}$ disjoint finite elements
$\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$ (of either quadrilateral
or triangular shape) in the following manner:

$\bigcup_{s=1}^{\prescript{\textrm{p}}{rg}{}N^{\textrm{c}}}\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}=\mathcal{H}_{{rg}},\quad\prescript{\textrm{p}}{rg}{}\Omega_{i}^{\textrm{c}}\cap\prescript{\textrm{p}}{rg}{}\Omega_{j}^{\textrm{c}}=\varnothing,~{}\forall
i\neq j$ (22)

Using the spatial mesh presented above, the nonlocal stress $\bm{\sigma}$ can
be expressed as:

$\bm{\sigma}(\bm{x}(\bm{\xi}_{rg}^{\textrm{p}}))=\bm{C}\left(\bm{x}(\bm{\xi}_{rg}^{\textrm{p}})\right):\sum_{s=1}^{\prescript{\textrm{p}}{rg}{}N^{\textrm{c}}}\int_{\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}}\mathcal{K}\left(\bm{x}(\bm{\xi}_{rg}^{\textrm{p}}),\bm{x}^{\prime}\right)\bm{\epsilon}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}$
(23)

Particular attention is required for mesh generation when considering weakly-
singular nonlocal kernel functions and Gauss-Jacobi quadrature. As discussed
in [40], integrals with weak singularity can be numerically evaluated by
Gauss-Jacobi quadrature only when the singular points coincide with the
boundary; see remarks in §2.5. In this regard, the singular point (for
example, $\bm{x}^{\prime}=\bm{x}$ in fractional-order nonlocal elasticity, see
§S1 of the SM) must be employed as a node in the child mesh such that it lies
on the boundary of the child elements
$\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$.

####  Numerical integration

In this section, we develop the approach to numerically evaluate the nonlocal
stress at the various quadrature points of the parent mesh. Similar to the
approach developed in §2.3.2, we first transform the child elements
$\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$ in 2D Cartesian
coordinates into 2D isoparametric elements
$\prescript{\textrm{p}}{rg}{}\hat{\Omega}_{s}^{\textrm{c}}$. By doing this, a
given point $\bm{x}=(x,y)$ in
$\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$, can now be mapped into:

$\bm{x}(\bm{\xi}_{s}^{\textrm{c}})=\sum_{d=1}^{d^{\textrm{c}}}\Phi_{d}^{\textrm{c}}(\bm{\xi}_{s}^{\textrm{c}})\prescript{\textrm{p}}{rg}{}\bm{x}_{sd}^{\textrm{c}}$
(24)

where $\Phi_{d}^{\textrm{c}}$ and $d^{\textrm{c}}$ are the shape function and
the number of nodes $\prescript{\textrm{p}}{rg}{}\bm{x}_{sd}^{\textrm{c}}$ in
child element $\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$,
respectively. Using the Gauss quadrature principles, the integral over
$\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}$ can be approximated in
isoparametric elements
$\prescript{\textrm{p}}{rg}{}\hat{\Omega}_{s}^{\textrm{c}}$, as:

$\int_{\prescript{\textrm{p}}{rg}{}\Omega_{s}^{\textrm{c}}}\mathcal{K}\left(\bm{x}(\bm{\xi}_{g}^{\textrm{p}}),\bm{x}^{\prime}\right)\bm{\epsilon}(\bm{x}^{\prime})\mathrm{d}\bm{x}^{\prime}=\sum_{h=1}^{g^{\textrm{c}}}w_{sh}^{\textrm{c}}|\prescript{\textrm{p}}{rg}{}\bm{J}_{s}^{\textrm{c}}|\mathcal{K}\left(\bm{x}(\bm{\xi}_{g}^{\textrm{p}}),\bm{x}^{\prime}(\bm{\xi}_{sh}^{\textrm{c}})\right)\bm{\epsilon}(\bm{\xi}_{sh}^{\textrm{c}})$
(25)

where $g^{\textrm{c}}$ is the number of Gauss points
$\bm{\xi}_{sh}^{\textrm{c}}$ and weight $w_{sh}^{\textrm{c}}$ in child
elements; $|\prescript{\textrm{p}}{rg}{}\bm{J}_{s}^{\textrm{c}}|$ is the
determinant of Jacobian matrix
$\prescript{\textrm{p}}{rg}{}\bm{J}_{s}^{\textrm{c}}$ defined by:

$\prescript{\textrm{p}}{rg}{}\bm{J}_{s}^{\textrm{c}}=\prescript{\textrm{p}}{rg}{}\frac{\partial\bm{x}}{\partial\bm{\xi}_{s}^{\textrm{c}}}=\sum_{d=1}^{d^{\textrm{c}}}\frac{\partial\Phi_{d}^{c}(\bm{\xi}_{s}^{\textrm{c}})}{\partial\bm{\xi}_{s}^{\textrm{c}}}\prescript{\textrm{p}}{rg}{}\bm{x}_{sd}^{\textrm{c}}=\sum_{d=1}^{d^{\textrm{c}}}\prescript{\textrm{p}}{rg}{}\bm{x}_{sd}^{\textrm{c}}\otimes\Phi_{d,\bm{\xi}_{s}^{\textrm{c}}}^{\textrm{p}}(\bm{\xi}_{s}^{\textrm{c}})$
(26)

##### Scale-Bridging: link to the parent domain.

Critical insights appear from a closer analysis of the child domain-based
nonlocal integral evaluation derived above. As evident from Eq. (25), the
determination of the stress requires an approximation of the local strain
$\bm{\epsilon}$ at the child Gauss points $\bm{\xi}_{sh}^{\textrm{c}}$. We
first discuss the physical relevance of this numerical expression and then
proceed to approximate the local strain. Note that, in the strain-driven
integral elasticity approach, the (infinitesimal) strain is assumed to be the
symmetric gradient of the displacement field, following local continuum
assumptions. The algorithmic translation of this definition, based on the fine
and coarse structural domains, indicates that the strain must be evaluated at
the parent domain. In fact, it is this connection that imparts an implicit
multiscale character to the overall model. Now, as the child mesh is decoupled
from the parent mesh, an added numerical mechanism is required for the
multiscale relay of response (here, strain).

We formulate a scale-bridging algorithm consisting of three different
transformations to numerically implement the multiscale relay of information,
described above. Without any loss of generality, we consider the approximation
of $\bm{\epsilon}$ at $\bm{\xi}_{sh}^{\textrm{c}}$ in the parent domain level.
From a high level perspective, for the scale-bridging of
$\bm{\xi}_{sh}^{\textrm{c}}$, we first find out the parent element
$\Omega_{\tilde{r}}^{\textrm{p}}$ ($\neq\Omega_{{r}}^{\textrm{p}}$, in
general; see details in step 2, below) that contains
$\bm{x}(\bm{\xi}_{sh}^{\textrm{c}})$, and then approximate
$\bm{\epsilon}(\bm{x}(\bm{\xi}_{sh}^{\textrm{c}}))$ within
$\Omega_{\tilde{r}}^{\textrm{p}}$. This scale-bridging requires the following
sequential transformations:

  * •

Step 1: to perform the inverse isoparametric transformation from child Gauss
point with isoparametric coordinates $\bm{\xi}_{sh}^{\textrm{c}}$ to the
Euclidean coordinates $\bm{x}(\bm{\xi}_{sh}^{\textrm{c}})$. Note that the
inverse isoparametric transformation can be achieved easily since we use
linear shape functions $\Phi_{d}^{\textrm{c}}$ for geometric interpolation.

  * •

Step 2: to determine the specific parent element
$\Omega_{\tilde{r}}^{\textrm{p}}$ that contains the child Gauss point
$\bm{x}(\bm{\xi}_{sh}^{\textrm{c}})$ (using its Cartesian coordinates). This
requires a search over all parent elements because the parent element
$\Omega_{\tilde{r}}^{\textrm{p}}$ containing $\bm{\xi}_{sh}^{\textrm{c}}$, is
more likely different from the reference parent element
$\Omega_{{r}}^{\textrm{p}}$ (that is used in this §2.4 to derive the child
domain approximation; see introduction to §2.4). In fact,
$\\{\Omega_{\tilde{r}}^{\textrm{p}},\Omega_{{r}}^{\textrm{p}}\\}$ can be
identified as a pair of interacting parent elements (or equivalently, the
discrete realization of the interacting pair $\\{\bm{x},\bm{x}^{\prime}\\}$ in
Eq. (1).

  * •

Step 3: to perform the forward isoparametric transformation (in the parent
domain) to obtain the isoparametric coordinates of the child Gauss point
$\bm{\xi}_{\tilde{r}}^{\textrm{p}}$ that belongs to the parent element
$\tilde{r}$.

The three-step scale-bridging, along with underlying mesh-decoupling, are
schematically illustrated in Fig. (1).

Figure 1: Schematic showing the finite element configuration including both
mesh-decoupling and scale-bridging algorithms. The overall discretization of
the square-shaped parent domain (or the material geometry) and the decoupled
child domain (or the nonlocal horizon). The parent domain is discretized
uniformly into parent meshes consisting of rectangular parent elements and
nodes (see the solid blue edges and blue nodes). The rectangular-shaped
nonlocal horizon for parent Gauss point $\bm{\xi}_{rg}^{\textrm{p}}$
(represented by the bold blue cross sign) in a given parent element
$\Omega_{r}^{\textrm{p}}$ (represented by the blue patch), is discretized into
individual rectangular child meshes with rectangular child elements and nodes
(see the solid red edges and red nodes). The child and parent meshes are fully
decoupled. The dimension of the nonlocal horizon in the four directions is
outlined by the green dashed curves. The child Gauss point
$\bm{\xi}_{sh}^{\textrm{c}}$ in the child element $\Omega_{s}^{\textrm{c}}$,
highlighted by the bold red cross sign, is located at the parent element
$\Omega_{\tilde{r}}^{\textrm{p}}$ (represented by the red patch). The three-
step scale-bridging algorithm that maps the child Gauss point
$\bm{\xi}_{sh}^{\textrm{c}}$ to the parent point
$\bm{\xi}_{\tilde{r}}^{\textrm{p}}$ is illustrated in the left, where child
and parent element are labeled as $\Omega_{s}^{\textrm{c}}$ and
$\Omega_{\tilde{r}}^{\textrm{p}}$, respectively. Similarly, the conversion
between Euclidean and isoparametric coordinates are tracked by $x$-$y$ and
$\xi$-$\eta$ in each element, respectively.

Following the above-outlined strategy,
$\bm{\epsilon}(\bm{x}(\bm{\xi}_{sh}^{\textrm{c}}))$ is approximated as:

$\bm{\epsilon}(\bm{x}(\bm{\xi}_{sh}^{\textrm{c}}))=\frac{1}{2}\left[\bm{\nabla}_{\bm{x}}\bm{u}(\bm{x}(\bm{\xi}_{\tilde{r}}^{\textrm{p}}))+\bm{\nabla}_{\bm{x}}^{\textrm{T}}\bm{u}(\bm{x}(\bm{\xi}_{\tilde{r}}^{\textrm{p}}))\right]$
(27)

Now, analogous to Eq. (16), the gradient of the displacement field at the
child Gauss point $\bm{\xi}_{sh}^{\textrm{c}}$ can be approximated within the
connected parent element $\bm{\xi}_{\tilde{r}}^{\textrm{p}}$, in the following
manner:

$\bm{\nabla}_{\bm{x}}\bm{u}(\bm{x}(\bm{\xi}_{\tilde{r}}^{\textrm{p}}))=\sum_{u=1}^{d^{\textrm{p}}}\bm{u}_{\tilde{r}u}\otimes\bm{\nabla}_{\bm{\xi}_{\tilde{r}}^{\textrm{p}}}\Psi_{u}^{\textrm{p}}(\bm{\xi}_{\tilde{r}}^{\textrm{p}})\frac{\partial\bm{\xi}_{\tilde{r}}^{\textrm{p}}}{\partial\bm{x}}=\sum_{u=1}^{d^{\textrm{p}}}\bm{u}_{\tilde{r}u}\otimes(\bm{J}_{\tilde{r}}^{\textrm{p}})^{-T}\bm{\nabla}_{\bm{\xi}_{\tilde{r}}^{\textrm{p}}}\Psi_{u}^{\textrm{p}}(\bm{\xi}_{\tilde{r}}^{\textrm{p}})$
(28)

where $\bm{u}_{\tilde{r}u}$ denotes DOF at nodal point
$u\in\\{1,2,...,d^{\mathrm{p}}\\}$ in the parent element
$\Omega_{\tilde{r}}^{\mathrm{p}}$. Note that, it is the inverse isoparametric
transformation that allows the mapping of the child Gauss point within the
child isoparametric coordinates and the Cartesian coordinates, such that the
displacement gradient at any child Gauss point can be approximated by using
DOF defined in parent elements.

###  Closing remarks on mesh-decoupling

The mesh-decoupling and scale-bridging algorithms deserve some additional
remarks. We emphasize that the mesh-decoupling essentially implies that the
child mesh (and the selection of its topology) is completely independent from
that of the parent mesh. Note that a physically-consistent multiscale transfer
of information is achieved via both the (three-step) scale-bridging approach,
and the theoretical thermodynamic consistency. The independence of the child
mesh topology allows the use of highly modular unstructured child domain
discretizations to account for different nonlocal kernels and irregular
horizons (as it will be shown in §5). The weakly-singular kernels and the
associated Gauss-Jacobi quadrature principle, for example, can be used only by
leveraging independent child mesh since the Gauss-Jacobi principle requires
the singularity point (that is the parent Gauss point
$\bm{\xi}_{rg}^{\textrm{p}}$) to be exactly at element edges, while in
comparison, if using parent mesh only (like [25]), the singular parent Gauss
point $\bm{\xi}_{rg}^{\textrm{p}}$ will never reach the edges in parent
meshes. We refer interested readers to §S1 of the supplementary material (SM)
for more discussion on the advantages of using the Gauss-Jacobi quadrature and
the independent child mesh.

Finally, a summary of the most important steps of $\text{M}^{2}$-FEM in
Algorithm 1 (with a specific focus on the mesh-decoupling technique), reveals
the scalable nature of the proposed concept. Evidently, the mesh-decoupling
technique can be extended in either directions (coarser or finer) from the
parent and/or child domains identified in Algorithm 1. The extension would
require an additional nested set of decoupled-meshes (analogous to that used
in Algorithm 1) and another scale-bridging to relay the relevant response
information. The resulting hierarchically stacked computational model can be
used to simulate highly multiscale problems. Although the need for such
hierarchically stacked multiscale theories has been highlighted (see [18, 24,
11]), we anticipate that the significant or prohibitive computational
resources that would be required to simulate these theories by existing finite
element methodologies, has prevented their development. Results presented in
§5 indicate that the proposed $\text{M}^{2}$-FEM could potentially enable the
development and implementation of hierarchically stacked multiscale theories.

Approximate weak-form by evaluating all variables defined at the coarse
(parent) scale

Approximate variables that are defined at fine (child) scale in weak-form

Mesh-decoupling for a given parent- and child-scale combination for  _every
element $\Omega_{r}^{\mathrm{p}}$ and corresponding quadrature points
$\bm{\xi}^{\mathrm{p}}_{r}$ in the parent mesh_ do

$\star$ Identify the child domain at the quadrature points
$\bm{\xi}^{\mathrm{p}}_{r}$ in $\Omega^{\mathrm{p}}_{r}$

$\star$ Generate an independent child mesh decoupled from parent mesh

$\star$ Numerically evaluate variables on child mesh using discretized
expressions

$\star$ Employ scale-bridging algorithm to relay information to parent domain

$\star$ Assemble all child domain approximations within weak-form integral

if  _variables defined at another (further) finer scale exist in the weak-
form_ then

$\quad-$ Update the child and the finer scales as the new parent and child
scales, respectively

$\quad-$ Repeat the parent-child mesh-decoupling approximation algorithm
above, recursively

Algorithm 1 Weak-form approximation of nested integrals in hierarchical
material scales via $\text{M}^{2}$-FEM

##  Implementation: matrix representation and assembled algebraic equations

The above two sections complete the development of the finite element
discretization and approximation in both the parent and child domains. Here
below, we derive the matrix representation of the different finite element
approximations and assemble them coherently to formulate the algebraic
equations of motion. It is common knowledge that these steps are essential
(and ideal) for an efficient computer implementation of the solver.

###  Background: matrix representations of primary variables and expressions

We consider the $r^{\textrm{th}}$ parent element $\Omega_{r}^{\textrm{p}}$,
the displacement DOF $\bm{u}$, and the weighting function $\delta\bm{v}$
defined on the element nodes. The finite element approximation of these
quantities can be expressed as the following vectors:

$\displaystyle\\{\bm{u}\\}_{r}^{\textrm{p}}$
$\displaystyle=\\{\\{\bm{u}_{x}\\}_{r}^{\mathrm{p}},\\{\bm{u}_{y}\\}_{r}^{\mathrm{p}}\\}^{\mathrm{T}}=\\{\\{\bm{u}_{x1},\bm{u}_{x2},...,\bm{u}_{xd^{\mathrm{p}}}\\},\\{\bm{u}_{y1},\bm{u}_{y2},...,\bm{u}_{yd^{\mathrm{p}}}\\}\\}^{\mathrm{T}}$
(29a) $\displaystyle\\{\delta\bm{v}\\}_{r}^{\mathrm{p}}$
$\displaystyle=\\{\\{\delta\bm{v}_{x}\\}_{r}^{\mathrm{p}},\\{\delta\bm{v}_{y}\\}_{r}^{\mathrm{p}}\\}^{\textrm{T}}=\\{\\{\delta\bm{v}_{x1},\delta\bm{v}_{x2},...,\delta\bm{v}_{xd^{\textrm{p}}}\\},\\{\delta\bm{v}_{y1},\delta\bm{v}_{y2},...,\delta\bm{v}_{yd^{\textrm{p}}}\\}\\}^{\textrm{T}}$
(29b)

in which $\bm{u}_{xr}$ and $\bm{u}_{yr}$ (or, $\delta\bm{v}_{xr}$ and
$\delta\bm{v}_{yr}$) are the previously defined solution vectors in size
$k\times 1$, that represent the nodal DOF (or, weighting functions) in the $x$
and $y$ direction, respectively. In that regard, the DOF
$\\{\bm{u}_{x}\\}_{r}^{\textrm{p}}$ and $\\{\bm{u}_{y}\\}_{r}^{\textrm{p}}$
(or, $\\{\delta\bm{v}_{x}\\}_{r}^{\textrm{p}}$ and
$\\{\delta\bm{v}_{y}\\}_{r}^{\textrm{p}}$) in a given parent element are all
row vectors of size $m^{\textrm{p}}\times 1$. Consequently, the length of
vectors $\\{\bm{u}\\}_{r}^{\textrm{p}}$ and
$\\{\delta\bm{v}\\}_{r}^{\textrm{p}}$ is $2m^{\textrm{p}}$. Similarly, the DOF
and weighting functions on the boundary element $\Gamma_{r}^{\textrm{t}}$,
both in size $2m^{\textrm{t}}\times 1$, can be expressed as:

$\displaystyle\\{\bm{u}\\}_{r}^{\textrm{t}}$
$\displaystyle=\\{\\{\bm{u}_{x}\\}_{r}^{\textrm{t}},\\{\bm{u}_{y}\\}_{r}^{\textrm{t}}\\}^{\textrm{T}}=\\{\\{\bm{u}_{x1},\bm{u}_{x2},...,\bm{u}_{xd^{\textrm{t}}}\\},\\{\bm{u}_{y1},\bm{u}_{y2},...,\bm{u}_{yd^{\textrm{t}}}\\}\\}^{\textrm{T}}$
(30a) $\displaystyle\\{\delta\bm{v}\\}_{r}^{\textrm{t}}$
$\displaystyle=\\{\\{\delta\bm{v}_{x}\\}_{r}^{\textrm{t}},\\{\delta\bm{v}_{y}\\}_{r}^{\textrm{t}}\\}^{\textrm{T}}=\\{\\{\delta\bm{v}_{x1},\delta\bm{v}_{x2},...,\delta\bm{v}_{xd^{\textrm{t}}}\\},\\{\delta\bm{v}_{y1},\delta\bm{v}_{y2},...,\delta\bm{v}_{yd^{\textrm{t}}}\\}\\}^{\textrm{T}}$
(30b)

We retrace the above steps for developing matrix representations of the shape
functions $\bm{\Phi}^{\textrm{p}}$ and $\bm{\Psi}^{\textrm{p}}$ in parent
elements, that perform the isogeometric transformation and the interpolation
of solution variables (also weighting functions, see Eq. (15)), respectively.
Following their definitions, we have:

$\displaystyle\\{\bm{\Phi}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\\}$
$\displaystyle=\\{\Phi_{1}(\bm{\xi}_{r}^{\textrm{p}}),\Phi_{2}(\bm{\xi}_{r}^{\textrm{p}}),...,\Phi_{d^{\textrm{p}}}(\bm{\xi}_{r}^{\textrm{p}})\\}^{\textrm{T}}$
(31a) $\displaystyle\\{\bm{\Psi}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\\}$
$\displaystyle=\\{\Psi_{1}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}}),\Psi_{2}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}}),...,\Psi_{m^{\textrm{p}}}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\\}^{\textrm{T}}$
(31b)

As we have clarified in §2.2, depending on the interpolation, the two sets of
shape functions can differ in both the definitions and dimensions (note that
$\\{\bm{\Phi}\\}^{\textrm{p}}$ is in size $d^{\textrm{p}}\times 1$ and
$\\{\bm{\Psi}^{\textrm{p}}\\}$ is in size $m^{\textrm{p}}\times 1$). In
analogy with the above formulation, shape functions $\bm{\Phi}^{\textrm{t}}$
and $\bm{\Psi}^{\textrm{t}}$ defined on boundary elements
$\Gamma_{r}^{\textrm{t}}$ can be also expressed as:

$\displaystyle\\{\bm{\Phi}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}})\\}$
$\displaystyle=\\{\Phi_{1}(\bm{\xi}_{r}^{\textrm{t}}),\Phi_{2}(\bm{\xi}_{r}^{\textrm{t}}),...,\Phi_{d^{\textrm{t}}}(\bm{\xi}_{r}^{\textrm{t}})\\}^{\textrm{T}}$
(32a) $\displaystyle\\{\bm{\Psi}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}})\\}$
$\displaystyle=\\{\Psi_{1}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}}),\Psi_{2}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}}),...,\Psi_{m^{\textrm{t}}}^{\textrm{t}}(\bm{\xi}_{r}^{\textrm{t}})\\}^{\textrm{T}}$
(32b)

We recast the original tensorial weak-form governing equations into the Voigt
matrix form. For nonlocal plane-strain problems with isotropic and homogeneous
elastic media, the second-order strain tensor, fourth-order stiffness tensor,
and the second-order nonlocal stress tensor can be recast into matrices as:

$[\bm{C}]_{3\times 3}=\begin{bmatrix}2\mu+\lambda&\lambda&0\\\
\lambda&2\mu+\lambda&0\\\ 0&0&2\mu\\\
\end{bmatrix},\quad\\{\bm{\sigma}\\}_{3\times 1}=\begin{bmatrix}\sigma_{xx}\\\
\sigma_{yy}\\\ \sigma_{xy}\end{bmatrix},\quad\\{\bm{\epsilon}\\}_{3\times
1}=\begin{bmatrix}\epsilon_{xx}\\\ \epsilon_{yy}\\\
\epsilon_{xy}\end{bmatrix}$ (33)

such that the nonlocal constitutive relation in Eq. (1) can be expressed as:

$\\{\bm{\sigma}\\}=[\bm{C}]\int_{\mathcal{H}}\mathcal{K}(\bm{x},\bm{x}^{\prime})\\{\bm{\epsilon}(\bm{x}^{\prime})\\}\mathrm{d}\bm{x}^{\prime}$
(34)

By using this technique, the double dot product
$\bm{\sigma}:\bm{\nabla}_{\bm{x}}\delta\bm{v}$ in the weak-form formulation
(see Eq. (4)) can be also recast into the following matrix representation:

$\bm{\sigma}:\bm{\nabla}_{\bm{x}}\delta\bm{v}=\\{\bm{\sigma}\\}^{\textrm{T}}\\{\bm{\nabla}_{\bm{x}}^{\textrm{sym}}\delta\bm{v}\\}_{3\times
1}=\\{\bm{\sigma}\\}^{\textrm{T}}\begin{bmatrix}\delta\bm{v}_{x,x}&\delta\bm{v}_{y,y}&\dfrac{1}{2}\left(\delta\bm{v}_{x,y}+\delta\bm{v}_{y,x}\right)\end{bmatrix}^{\textrm{T}}$
(35)

where $\bm{\nabla}_{\bm{x}}^{\textrm{sym}}$ denotes the symmetric gradient in
Voigt notation. By using the above framework, we recast the $\text{M}^{2}$-FEM
algorithm into an algebraic representation that can be easily translated into
a scientific computing program.

###  Matrix representation of numerically evaluated approximations

We first evaluate the terms containing the weighting functions $\delta\bm{v}$.
Specifically, we start with the finite element approximation of the weighting
functions $\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})$ at Gauss points
$\bm{\xi}_{rg}^{\textrm{p}}$. By combining the numerical scheme in Eq. (15)
and the vectorial form of shape functions $\bm{\Psi}^{\textrm{p}}$, we obtain
the matrix representation of $\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})$ as
follows:

$\\{\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})\\}=\\{\bm{L}\\}_{1\times
2m^{\textrm{p}}}^{(\bm{\Psi}^{\textrm{p}})}\\{\delta\bm{v}\\}_{r}=\\{\bm{\Psi}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}}),\bm{\Psi}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\\}^{\textrm{T}}\\{\delta\bm{v}\\}_{r}$
(36)

Note that we use $\\{\bm{\Psi}^{\textrm{p}}(\bm{\xi}_{r}^{\textrm{p}})\\}$
twice in the coefficient matrix $\\{\bm{L}\\}^{(\bm{\Psi}^{\textrm{p}})}$ to
account for the approximation of both ${\delta\bm{v}_{x}}_{r}$ and
${\delta\bm{v}_{y}}_{r}$. The other weighting function related term in the
parent element, which is the gradient of weighting function
$\bm{\nabla}_{\bm{x}}\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})$ in Eq. (16),
can be replaced by its vectorial form
$\\{\bm{\nabla}_{\bm{x}}^{\textrm{sym}}\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})\\}$:

$\\{\bm{\nabla}_{\bm{x}}^{\textrm{sym}}\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{p}})\\}=[\bm{G}]_{3\times
2m^{\textrm{p}}}^{(\delta\bm{v})}\\{\delta\bm{v}\\}_{r}=\frac{1}{2}\begin{bmatrix}\\{\bm{G}_{rx}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{\textrm{T}}&\\{\bm{0}\\}_{1\times
m^{\textrm{p}}}\\\\[4.0pt] \\{\bm{0}\\}_{1\times
m^{\textrm{p}}}&\\{\bm{G}_{ry}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{\textrm{T}}\\\\[4.0pt]
\\{\bm{G}_{rx}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{\textrm{T}}&\\{\bm{G}_{ry}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{\textrm{T}}\end{bmatrix}\\{\delta\bm{v}\\}_{r}$
(37)

where $[\bm{G}]_{3\times 2m^{\textrm{p}}}^{(\delta\bm{v})}$ is the coefficient
matrix in size $3\times 2m^{\textrm{p}}$; $\\{\bm{0}\\}_{1\times
m^{\textrm{p}}}$ is the zero vector in size $1\times m^{\textrm{p}}$;
$\\{\bm{G}_{rx}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}$ and
$\\{\bm{G}_{rx}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}$ are two
coefficient vectors in size $m^{\textrm{p}}\times 1$ that represent the first-
order spatial gradient along $x$ and $y$ directions. The two gradient
coefficient vectors can be evaluated by using the matrix
$[\bm{G}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})]$ (see Eq. (16)):

$[\\{\bm{G}_{rx}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\},\\{\bm{G}_{ry}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}]=[\bm{G}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})]_{m^{\textrm{p}}\times
2}=[\bm{\nabla}_{\bm{\xi}_{r}^{\textrm{p}}}\Psi_{v}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})(\bm{J}_{rg}^{\textrm{p}})^{-1}]=\begin{bmatrix}\bm{\nabla}_{\bm{\xi}_{r}^{\textrm{p}}}\Psi_{1}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})(\bm{J}_{rg}^{\textrm{p}})^{-1}\\\\[4.0pt]
\bm{\nabla}_{\bm{\xi}_{r}^{\textrm{p}}}\Psi_{2}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})(\bm{J}_{rg}^{\textrm{p}})^{-1}\\\\[4.0pt]
\vdots\\\\[4.0pt]
\bm{\nabla}_{\bm{\xi}_{r}^{\textrm{p}}}\Psi_{m^{\textrm{p}}}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})(\bm{J}_{rg}^{\textrm{p}})^{-1}\end{bmatrix}$
(38)

Finally, the matrix form of the approximated $\delta\bm{v}$ on boundary
elements $\Gamma_{r}^{\textrm{t}}$ can be readily given by:

$\\{\delta\bm{v}(\bm{\xi}_{rg}^{\textrm{t}})\\}=\\{\bm{L}(\bm{\xi}_{rg}^{\textrm{t}})\\}_{1\times
2m^{\textrm{t}}}^{(\bm{\Psi}^{\textrm{t}})}\\{\delta\bm{v}\\}_{r}^{\textrm{t}}=\\{\bm{\Psi}^{\textrm{t}}(\bm{\xi}_{rg}^{\textrm{t}}),\bm{\Psi}^{\textrm{t}}(\bm{\xi}_{rg}^{\textrm{t}})\\}^{\textrm{T}}\\{\delta\bm{v}\\}_{r}^{\textrm{t}}$
(39)

We then derive the matrix representation for terms containing $\bm{u}$.
Similar to Eq. (36), the matrix representation of
$\ddot{\bm{u}}(\bm{\xi}_{rg}^{\textrm{p}})$ shown in Eq. (15) can be obtained
immediately as:

$\\{\ddot{\bm{u}}(\bm{\xi}_{rg}^{\textrm{p}})\\}=\\{\bm{L}(\bm{\xi}_{rg}^{\textrm{p}})\\}_{1\times
2m^{\textrm{p}}}^{(\bm{\Psi}^{\textrm{p}})}\\{\ddot{\bm{u}}\\}_{r}$ (40)

The matrix formulation with respect to the other solution variable related
term $\bm{\sigma}$ is more complicated due to the nonlocal integration. To
formulate the nonlocal stress $\bm{\sigma}$, we start with the local strain
tensor $\bm{\epsilon}(\bm{x}(\bm{\xi}_{sh}^{\textrm{c}}))$ defined at the
child Gauss point $\bm{\xi}_{sh}^{\textrm{c}}$. In analogy with
$\\{\bm{\nabla}_{\bm{x}}^{\textrm{sym}}\delta\bm{v}\\}$ in Eq. (37), the
second-order strain tensor $\bm{\epsilon}$ can be also reformulated
accordingly into the following Voigt form:

$\\{\bm{\epsilon}(\bm{x}(\bm{\xi}_{sh}^{\textrm{c}}))\\}=[\bm{G}(\bm{\xi}_{sh}^{\textrm{c}})]_{3\times
2m^{\textrm{p}}}^{(\bm{u})}\\{\bm{u}\\}_{r}=\frac{1}{2}\begin{bmatrix}\\{\bm{G}_{\tilde{r}x}^{\textrm{p}}(\bm{\xi}_{\tilde{r}g}^{\textrm{p}})\\}^{\textrm{T}}&\\{\bm{0}\\}_{1\times
m^{\textrm{p}}}\\\\[4.0pt] \\{\bm{0}\\}_{1\times
m^{\textrm{p}}}&\\{\bm{G}_{\tilde{r}y}^{\textrm{p}}(\bm{\xi}_{\tilde{r}g}^{\textrm{p}})\\}^{\textrm{T}}\\\\[4.0pt]
\\{\bm{G}_{\tilde{r}y}^{\textrm{p}}(\bm{\xi}_{\tilde{r}g}^{\textrm{p}})\\}^{\textrm{T}}&\\{\bm{G}_{\tilde{r}x}^{\textrm{p}}(\bm{\xi}_{\tilde{r}g}^{\textrm{p}})\\}^{\textrm{T}}\end{bmatrix}\\{\bm{u}\\}_{\tilde{r}}$
(41)

in which the
$\\{\bm{G}_{\tilde{r}x}^{\textrm{p}}(\bm{\xi}_{\tilde{r}g}^{\textrm{p}})\\}$
and
$\\{\bm{G}_{\tilde{r}y}^{\textrm{p}}(\bm{\xi}_{\tilde{r}g}^{\textrm{p}})\\}$
share the same definition with their counterparts, namely the
$\\{\bm{G}_{rx}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}$ and
$\\{\bm{G}_{rx}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})\\}$, that are defined
in Eq. (38). Note that $\bm{\xi}_{sh}^{\textrm{c}}$ and
$\bm{\xi}_{\tilde{r}g}^{\textrm{p}}$ are essentially the same child Gauss
point, but represented in the child and parent mesh, respectively (see the
discussion in §2.4.2). Further substituting Eq. (41) into the nonlocal
integral, we obtain the matrix representation of the nonlocal stress
$\bm{\sigma}$:

$\bm{\sigma}(\bm{x}(\bm{\xi}_{rg}^{\textrm{p}}))=[\bm{E}(\bm{\xi}_{rg}^{\textrm{p}})]_{3\times
2kN^{\textrm{s}}}\\{\bm{u}\\}_{2kN^{\textrm{s}}\times 1}$ (42)

with the coefficient matrix $[\bm{E}(\bm{\xi}_{rg}^{\textrm{p}})]$ in size
$3\times 2kN^{\textrm{s}}$ is calculated by:

$[\bm{E}(\bm{\xi}_{rg}^{\textrm{p}})]=[\bm{C}(\bm{x}(\bm{\xi}_{rg}^{\textrm{p}}))]\left[\sum_{s=1}^{\prescript{\textrm{p}}{rg}{}N^{\textrm{c}}}\sum_{h=1}^{n^{\textrm{c}}}\omega_{sh}^{\textrm{c}}|\prescript{\textrm{p}}{rg}{}\bm{J}_{s}^{\textrm{c}}|\mathcal{K}\left(\bm{x}(\bm{\xi}_{g}^{\textrm{p}}),\bm{x}^{\prime}(\bm{\xi}_{sh}^{\textrm{c}})\right)[\bm{G}(\bm{\xi}_{sh}^{\textrm{c}})]^{(\bm{u})}\right]$
(43)

Note that, due to the nonlocal integration, the coefficient matrix
$[\bm{E}(\bm{\xi}_{rg}^{\textrm{p}})]$ for the nonlocal stress $\bm{\sigma}$
is not in shape $3\times 2m^{\textrm{p}}$ any more; instead, it requires
displacement information from all the nodes within the nonlocal horizon at the
current parent Gauss point $\bm{\xi}_{rg}^{\textrm{p}}$. Depending on the size
of nonlocal horizon, the total number of nodes interacting with the current
parent Gauss point varies. For this purpose, in Eq. (42) we define the
coefficient matrix to be in size $3\times 2kN^{\textrm{s}}$ such that it
allows including potential nonlocal interactions from all the $N^{\textrm{s}}$
nodes in the parent domain. Accordingly, we use the global DOF vector
$\\{\bm{u}\\}$ (which is in size $2kN^{\textrm{s}}\times 1$) to denote DOF at
all the $N^{\textrm{s}}$ nodes.

###  Assembly and derivation of algebraic equations of motion

The matrix representation of the element-wise weak-form integration now can be
derived readily. By assembling all matrices into Eqs. (17,21), the weak-form
integration can be reorganized into the following matrix forms:

$\displaystyle\int_{\Omega_{r}^{\textrm{p}}}\bm{\sigma}:\bm{\nabla}_{\bm{x}}\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}$
$\displaystyle=\\{\bm{u}\\}^{\textrm{T}}[\bm{K}_{r}^{\textrm{p}}]_{2kN^{\textrm{s}}\times
2m^{\textrm{p}}}\\{\delta\bm{v}\\}_{r}^{\textrm{p}}$ (44a)
$\displaystyle\int_{\Omega_{r}^{\textrm{p}}}\rho(\bm{f}-\ddot{\bm{u}})\delta\bm{v}\mathrm{d}\Omega_{r}^{\textrm{p}}$
$\displaystyle=\\{\bm{F}_{r}^{\textrm{p}}\\}_{1\times
2m^{\textrm{p}}}\\{\delta\bm{v}\\}_{r}^{\textrm{p}}-\left[\\{\ddot{\bm{u}}\\}_{r}^{\textrm{p}}\right]^{\textrm{T}}[\bm{M}_{r}^{\textrm{p}}]_{2m^{\textrm{p}}\times
2m^{\textrm{p}}}\\{\delta\bm{v}\\}_{r}^{\textrm{p}}$ (44b)
$\displaystyle\int_{\Gamma_{r}^{\textrm{t}}}\overline{\bm{t}}\delta\bm{v}\mathrm{d}\Gamma_{r}^{\textrm{t}}$
$\displaystyle=\\{\bm{F}_{r}^{\textrm{t}}\\}_{1\times
2m^{\textrm{t}}}\\{\delta\bm{v}\\}_{r}^{\textrm{t}}$ (44c)

where $[\bm{K}_{r}^{\textrm{p}}]$, $[\bm{M}_{r}^{\textrm{p}}]$,
$\\{\bm{F}_{r}^{\textrm{p}}\\}$, and $\\{\bm{F}_{r}^{\textrm{t}}\\}$ are
coefficient matrices and vectors that are defined, respectively, as:

$\displaystyle[\bm{K}_{r}^{\textrm{p}}]$
$\displaystyle=\sum_{g=1}^{g^{\textrm{p}}}{w}_{rg}^{\textrm{p}}|\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|[\bm{E}(\bm{\xi}_{rg}^{\textrm{p}})]^{\textrm{T}}[\bm{G}(\bm{\xi}_{rg}^{\textrm{p}})]^{(\delta\bm{v})}$
(45a) $\displaystyle[\bm{M}_{r}^{\textrm{p}}]$
$\displaystyle=\sum_{g=1}^{g^{\textrm{p}}}w_{rg}^{\textrm{p}}|\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|\rho(\bm{\xi}_{rg}^{\textrm{p}})\left[\\{\bm{L}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{(\bm{\Psi}^{\textrm{p}})}\right]^{\textrm{T}}\\{\bm{L}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{(\bm{\Psi}^{\textrm{p}})}$
(45b) $\displaystyle\\{\bm{F}_{r}^{\textrm{p}}\\}$
$\displaystyle=\sum_{g=1}^{g^{\textrm{p}}}w_{rg}^{\textrm{p}}|\bm{J}_{r}^{\textrm{p}}(\bm{\xi}_{rg}^{\textrm{p}})|\rho(\bm{\xi}_{rg}^{\textrm{p}})\bm{f}(\bm{\xi}_{rg}^{\textrm{p}})\\{\bm{L}(\bm{\xi}_{rg}^{\textrm{p}})\\}^{(\bm{\Psi}^{\textrm{p}})}$
(45c) $\displaystyle\\{\bm{F}_{r}^{\textrm{t}}\\}$
$\displaystyle=\sum_{g=1}^{g^{\textrm{t}}}w_{rg}^{\textrm{t}}|\bm{J}_{r}^{\textrm{t}}((\bm{\xi}_{rg}^{\textrm{t}}))|\overline{\bm{t}}(\bm{\xi}_{rg}^{\textrm{t}})\\{\bm{L}(\bm{\xi}_{rg}^{\textrm{t}})\\}^{(\bm{\Psi}^{\textrm{t}})}$
(45d)

Iterating the above process for integration over all the parent and boundary
elements, we obtain the final matrix representation of the finite element
approximation to the weak-form governing equations in Eq. (4):

$\\{\ddot{\bm{u}}\\}^{\textrm{T}}[\bm{M}]_{2kN^{\textrm{s}}\times
2kN^{\textrm{s}}}\\{\delta\bm{v}\\}+\\{\bm{u}\\}^{\textrm{T}}[\bm{K}]_{2kN^{\textrm{s}}\times
2kN^{\textrm{s}}}\\{\delta\bm{v}\\}=\\{\bm{F}\\}_{1\times
2kN^{\textrm{s}}}\\{\delta\bm{v}\\}$ (46)

by assembling all the coefficient matrices and vectors
$[\bm{K}_{r}^{\textrm{p}}]$, $[\bm{M}_{r}^{\textrm{p}}]$,
$\\{\bm{F}_{r}^{\textrm{p}}\\}$, and $\\{\bm{F}_{r}^{\textrm{t}}\\}$ in each
elements together. To guarantee that Eq. (46) holds under any variation of
weighting function $\\{\delta\bm{v}\\}$, we must have:

$[\bm{M}]^{\textrm{T}}\\{\ddot{\bm{u}}\\}+[\bm{K}]^{\textrm{T}}\\{\bm{u}\\}=\\{\bm{F}\\}^{\textrm{T}}$
(47)

where $[\bm{M}]^{\textrm{T}}$, $[\bm{K}]^{\textrm{T}}$, and
$\\{\bm{F}\\}^{\textrm{T}}$ are the final mass matrix, the stiffness matrix,
and the force vector, respectively (note that the transpose operators are
applied for convenience). Equation (47) is the finally approximated algebraic
linear system equivalent of the integral nonlocal elasticity problem presented
in Eq. (10).

##  Validation and convergence studies

In this section, we investigate the validity and convergence of the developed
algorithm by simulating a set of 2D nonlocal plane strain problems. A common
parent domain, with size $L_{x}\times L_{y}=L\times L=1\textrm{m}\times
1\textrm{m}$ and Lamé’s parameter $\mu=1\textrm{Pa}$ and
$\lambda=1\textrm{Pa}$, is adopted for all the simulations. The material is
chosen to be isotropic in nature following the assumptions of Eringen’s theory
[17, 13]. In order to investigate the performance of $\text{M}^{2}$-FEM, we
define and perform two parametric case studies by varying the definition of
the kernel and the size of the horizon.

###  Definitions of case studies

Two case studies are adopted:

  * •

Case study 1 - Non-singular bi-exponential kernel adopted in [17] is assumed:

$\mathcal{K}(x,y,x^{\prime},y^{\prime},\tau)=\dfrac{1}{\pi\tau}\mathrm{exp}\left(-\dfrac{|x-x^{\prime}|^{2}}{\tau}\right)\mathrm{exp}\left(-\dfrac{|y-y^{\prime}|^{2}}{\tau}\right),\quad(x^{\prime},y^{\prime})\in\mathcal{H}(x,y)$
(48)

where $\tau$ is the nonlocal strength parameter that controls the decay of the
bi-exponential nonlocal kernel function.

  * •

Case study 2 - Weakly-singular power-law kernel adopted in fractional-order
nonlocal theories [11, 12], is assumed:

$\mathcal{K}(x,y,x^{\prime},y^{\prime},\alpha)=\frac{1}{\Gamma^{2}(1-\alpha)}\frac{1}{|x-x^{\prime}|^{\alpha}|y-y^{\prime}|^{\alpha}},\quad(x^{\prime},y^{\prime})\in\mathcal{H}(x,y)$
(49)

where $\alpha\in(0,1)$ controls the spatial decay of the kernel and is called
as ’fractional-order’, and $\Gamma(\cdot)$ denotes the Gamma function [30].
Note that the power-law kernel is (weakly) singular along lines $x^{\prime}=x$
and $y^{\prime}=y$.

Figure (S1) in the SM illustrates the difference between the above kernel
functions. Note that, both the kernel functions defined above are positive-
definite and symmetric such that the nonlocal model is well-posed (see §2.1 or
[16]). In each case study, the nonlocal horizon is assumed to be isotropic and
symmetric in shape such that, at a given point $(x,y)$ in the parent domain,
the 2D nonlocal horizon can be expressed as:

$\mathcal{H}(x,y)=[\textrm{min}(x-l^{\textrm{nl}},0),\textrm{max}(x+l^{\textrm{nl}},L)]\times[\textrm{min}(y-l^{\textrm{nl}},0),\textrm{max}(y+l^{\textrm{nl}},L)]$
(50)

where $l^{\textrm{nl}}$ denotes the length-scale in the child domain. Note
that the minimum and maximum operators are required to contain the nonlocal
horizon within the material boundaries.

For the present analyses, the parent and child domains are assumed to be
squares for two specific reasons. First, this choice enables the derivation of
closed-form analytical expressions for specific loading conditions which are
utilized for validation. Second, this choice allows the definition of commonly
used computational resolution measures which provide clear insights on the
convergence characteristics of the integral solver. Next, the parent domain is
uniformly discretized with elements of size
$\Delta_{x}^{\textrm{p}}=\Delta_{y}^{\textrm{p}}=\Delta^{\textrm{p}}$, such
that the parent mesh number is
$M_{x}^{\textrm{p}}=M_{y}^{\textrm{p}}=M^{\textrm{p}}=L/\Delta^{\textrm{p}}$.
Unlike the parent mesh, the child mesh is unstructured due to the truncation
of nonlocal horizons at boundaries, and also the possible assumption of
weakly-singular nonlocal kernels (recall from §3.2, in this case the Gauss
quadrature points must coincide with the nodes). We have discussed the child
domain discretization strategy in §1 of the SM, in the interest of keeping the
manuscript concise. Given the non-uniform child mesh, we use the spatial
average of the child element size, denoted as $\Delta^{\textrm{c}}$, to
represent the child element size. The ratio of the parent and child element
sizes is denoted as
$r^{\textrm{p-c}}(=\Delta^{\textrm{p}}/\Delta^{\textrm{c}})$, and called as
relative resolution.

###  Validation studies

We validate the $\text{M}^{2}$-FEM by comparing the results of the solver
against the exact solution (the analytical expressions of the displacement
field) of the defined set (§4.1) of 2D plane strain integral elasticity
problems with spatially-varying loads. Since, analytical solutions of differ-
integral partial differential equations are not derived readily (even for very
simple loading conditions), we adopt a different strategy, as discussed in the
following. First, we assume the 2D displacement field to obey the following
analytical variation:

$\bm{u}(x,y)\equiv\\{u_{x}(x,y),u_{y}(x,y)\\}=\\{(-x^{2}+Lx)(-y^{2}+Ly),0\\}$
(51)

The assumed field satisfies homogeneous displacement boundary conditions, and
the bi-parabolic $u_{x}$ displacement achieves its maximum value
$\bar{u}_{x}^{\textrm{ref}}=0.0625$[m] at $(L/2,L/2)$. The above reference
displacement field is substituted within the nonlocal model in §2.1 to obtain
the analytical expressions of the body loads ($f_{x}$ and $f_{y}$); see SM §2.
Finally, the derived loads are used as input to $\text{M}^{2}$-FEM to simulate
the case studies and evaluate the accuracy of the numerical results (against
the reference solution in Eq. (51)).

Using the above strategy, the error noted in the simulation of a set of
parametric configurations, obtained by varying the nonlocal horizon
($l^{\textrm{nl}}$) and the kernel decay-rate parameter ($\tau$ for case study
and $\alpha$ for case study 2), are evaluated. Additionally, for a
comprehensive analysis, the resolution of either (a) parent mesh or (b) child
mesh is varied. The results are are presented in Table 1 and Table 2 (for the
case study 1 and case study 2, respectively) in terms of the simulated maximum
displacement (in the $x$-direction) and the relative error in terms of the
absolute difference between the simulated value from the reference value
($\bar{u}_{x}^{\textrm{ref}}=0.0625$[m]). In analyzing the impact of the child
mesh resolution, the parent resolution is fixed and the relative resolution
($r^{\textrm{p-c}}$) is varied. Note that the relative resolution is assumed
$r^{\textrm{p-c}}\sim\mathcal{O}(1)$ to prevent erroneous conclusions
resulting either due to computational instabilities (when $r^{\textrm{p-c}}\ll
1$) or due to poor efficiency (when $r^{\textrm{p-c}}\gg 1$). The full-field
simulation for a specific combination of the two kernels, nonlocal horizon
length, and mesh sizes, are presented in Fig. (2) and Fig. (3). The excellent
match between the simulated and reference responses successfully validates the
implementation of M2-FEM (particularly, the mesh-decoupling and scale-bridging
schemes, that are the pillars of M2-FEM).

Table 1: Relative error in the maximum displacement along $x$-direction obtained by using $\text{M}^{2}$-FEM, against the reference value, parameterized for different values of the nonlocal horizon length ($l^{nl}$) and kernel decay rate ($\tau$) for Case study 1. The effect of the (a) parent domain resolution and (b) relative resolution on accuracy are also analyzed. $l^{\textrm{nl}}$ | (a) $M^{\textrm{p}}$ | $\tau=0.002L$ | $\tau=0.003L$ | (b) $r^{\textrm{p-c}}$ | $\tau=0.002L$ | $\tau=0.003L$
---|---|---|---|---|---|---
$\overline{u}_{x}\times 100$ | Er[%] | $\overline{u}_{x}\times 100$ | Er[%] | $\overline{u}_{x}\times 100$ | Er[%] | $\overline{u}_{x}\times 100$ | Er[%]
$l^{\textrm{nl}}=0.2L$ | $8$ | 5.96 | 4.67 | 6.0647 | 2.96 | $1$ | 6.23 | 0.40 | 6.24 | 0.24
$16$ | 6.23 | 0.28 | 6.2450 | 0.08 | $2$ | 6.24 | 0.17 | 6.24 | 0.23
$24$ | 6.24 | 0.11 | 6.2377 | 0.12 | $3$ | 6.24 | 0.17 | 6.24 | 0.23
$l^{\textrm{nl}}=0.3L$ | $8$ | 5.96 | 4.65 | 6.087 | 2.64 | $1$ | 6.23 | 0.40 | 6.24 | 0.24
$16$ | 6.23 | 0.28 | 6.24 | 0.08 | $2$ | 6.24 | 0.17 | 6.24 | 0.23
$24$ | 6.24 | 0.11 | 6.4 | 0.20 | $3$ | 6.24 | 0.17 | 6.24 | 0.23
Table 2: Relative error in the maximum displacement along $x$-direction obtained by using $\text{M}^{2}$-FEM, against the reference value, parameterized for different values of the nonlocal horizon length ($l^{\textrm{nl}}$) and power-law orders ($\alpha$) for Case study 2. The effect of the (a) parent domain resolution and (b) relative resolution on accuracy are also analyzed. $l^{\textrm{nl}}$ | (a) $M^{\textrm{p}}$ | $\alpha=0.6$ | $\alpha=0.3$ | (b) $r^{\textrm{p-c}}$ | $\alpha=0.6$ | $\alpha=0.3$
---|---|---|---|---|---|---
$\overline{u}_{x}\times 100$ | Er[%] | $\overline{u}_{x}\times 100$ | Er[%] | $\overline{u}_{x}\times 100$ | Er[%] | $\overline{u}_{x}\times 100$ | Er[%]
$l^{\textrm{nl}}=0.3L$ | $8$ | 6.15 | 1.52 | 6.24 | 0.15 | $1$ | 6.26 | 0.19 | 6.22 | 0.43
$16$ | 6.20 | 0.75 | 6.28 | 0.45 | $2$ | 6.22 | 0.34 | 6.21 | 0.68
$24$ | 6.30 | 0.80 | 6.13 | 1.87 | $3$ | 6.223 | 0.43 | 6.21 | 0.70
$l^{\textrm{nl}}=0.5L$ | $8$ | 6.34 | 1.51 | 6.48 | 3.60 | $1$ | 6.28 | 0.42 | 6.35 | 1.64
$16$ | 6.27 | 0.40 | 6.32 | 1.19 | $2$ | 6.23 | 0.24 | 6.27 | 0.34
$24$ | 6.26 | 0.18 | 6.29 | 0.59 | $3$ | 6.23 | 0.36 | 6.25 | 0.08
Figure 2: $\text{M}^{2}$-FEM simulation of the nonlocal elastic configuration
with non-singular bi-exponential kernel in Case study 1. The nonlocal response
is presented in terms of the displacement field along the a $x-$direction
($u_{y}$) and b $x-$direction ($u_{y}$). The nonlocal model parameters are set
to $\tau_{x}=\tau_{y}=0.002$ and $l_{x}=l_{y}=0.5$, and the mesh parameters
are set to $\Delta_{x}^{\textrm{p}}=\Delta_{y}^{\textrm{p}}=1/30$ and
$\Delta_{x}^{\textrm{c}}=\Delta_{y}^{\textrm{c}}=1/30$. Both the weak-form
integral in the parent domain and the bi-exponential nonlocal integral in the
child domain are evaluated numerically using Gauss-Legendre quadrature. The
second-order Lagrange element is used for interpolation. Figure 3:
$\text{M}^{2}$-FEM simulation of the nonlocal elastic configuration with
singular (fractional-order) power-law kernel in Case study 2. The nonlocal
response is presented in terms of the displacement field along the a
$x-$direction ($u_{y}$) and b $x-$direction ($u_{y}$). The nonlocal model
parameters are set as $\alpha_{x}=\alpha_{y}=0.5$ and $l_{x}=l_{y}=0.5$, and
the mesh parameters are set as
$\Delta_{x}^{\textrm{p}}=\Delta_{y}^{\textrm{p}}=1/30$ and
$\Delta_{x}^{\textrm{c}}=\Delta_{y}^{\textrm{c}}=1/30$. The weak-form integral
is numerically evaluated by the Gauss-Legendre quadrature and the weakly-
singular nonlocal integral is evaluated by using Gauss-Jacobi quadrature. The
first-order Lagrange element is for interpolation.

###  Convergence analysis

Analogous to existing integral solvers [25, 26, 27, 28], the convergence of
the $\text{M}^{2}$-FEM involves a more elaborate analysis (when compared to
classical FEM) due to the additional child domain integral approximation that
is nested within the weak-form integral. In addition to completeness and
compatibility conditions typical of the finite element simulation of local
elasticity problems [39], the convergence of the nested parent-child integrals
in the simulation of integral elasticity theories is affected also by the
relative resolution of the parent and child meshes (equivalent to the dynamic
rate of convergence in mesh-coupled solvers [26, 27]; shown below) and the
nature of the kernel function (more particularly, the spatial decay rate and
its resolution within the child mesh) [14]. We present comprehensive
parametric studies by leveraging both h-refinement and the p-refinement to
investigate the convergence conditions in terms of the aforementioned factors.
In the subsequent analysis, we define the convergence criterion as a
(arbitrary and commonly adopted) threshold of 2% relative difference in the
Euclidean distance between the displacement fields simulated from successive
spatial refinements (in either the parent or child domain); the relative
difference is denoted as $\Delta$ simply, without any super- or sub-scripts
(to avoid confusion with the parent and child element sizes).

#####  h-refinement

We perform h-refinement in the child domain because the analysis of accuracy
presented in §4.2, as well as convergence studies conducted for existing
integral solvers [27, 28], indicate that the performance of the integral
solver is sensitive to the child domain resolution. In fact, we first show
that the 2D child domain resolution is equivalent to the dynamic rate of
convergence that explicitly determines the convergence in mesh-coupled
solvers:

$\left[\frac{l^{\mathrm{nl}}}{\Delta^{\mathrm{c}}_{x}}\right]\times\left[\frac{l^{\mathrm{nl}}}{\Delta^{\mathrm{c}}_{y}}\right]\equiv\frac{1}{r^{\textrm{p-c}}}\left[\frac{l^{\mathrm{nl}}}{\Delta^{\mathrm{p}}_{x}}\right]\times\frac{1}{r^{\textrm{p-c}}}\left[\frac{l^{\mathrm{nl}}}{\Delta^{\mathrm{p}}_{y}}\right]\xRightarrow{r^{\textrm{p-c}}=1}\underbrace{\left[\frac{l^{\mathrm{nl}}}{\Delta^{\mathrm{p}}_{x}}\right]}_{\mathcal{N}_{x}}\times\underbrace{\left[\frac{l^{\mathrm{nl}}}{\Delta^{\mathrm{p}}_{y}}\right]}_{\mathcal{N}_{y}}$
(52)

In the second equivalence above, we assume $r^{\textrm{p-c}}=1$, as is the
case for mesh-coupled solvers, and $\mathcal{N}_{x}\times\mathcal{N}_{y}$ is
the dynamic rate of convergence. Remarkably, this mathematical deduction
suggests that the overall convergence (in addition to the desired accuracy)
can be directly achieved by independently tuning the child mesh resolution
while still maintaining a coarse parent mesh for superior efficiency (thanks,
to the decoupling; see §5.2).

The results of a parametric convergence analysis corresponding to the case
studies introduced previously is presented in Table 3. As evident from the
results, the convergence for the exponential kernel (case study 1) is achieved
at $\mathcal{N}_{x}\times\mathcal{N}_{y}=6\times 6$ and for the power-law
kernel at $\mathcal{N}_{x}\times\mathcal{N}_{y}=5\times 5$ (case study 2). In
other terms, atleast 6 and 5 elements are required (on an average) along each
direction to achieve a detailed resolution of the displacement profile
following the exponential and power-law kernels, respectively. The (slightly)
slower convergence observed in the exponential kernel, when compared to the
power-law kernel, on refinement is due to the faster spatial decay of the
exponential kernel which requires a more detailed mesh for better resolution;
see similar analysis also in [14]. However, in comparison to existing solvers
which attain convergence around $\mathcal{N}_{\square}\sim\mathcal{O}(10)$
(where, $\square\in\\{x,y\\}$) [28, 25, 27, 26], $\text{M}^{2}$-FEM rapid
attains convergence on successive mesh refinement. The latter remarkable
outcome is a direct impact of the powerful ability of the decoupled child mesh
in enforcing a convergence (from the top, as in local FEM [39]) to the real
solution independently of the parent mesh; due to the minimal transmission of
approximation error to and from the parent mesh.

Table 3: Convergence of the $\text{M}^{2}$-FEM simulation of the nonlocal
problems in the two case studies via h-refinement.

Case study 1 | Case study 2
---|---
$l^{\textrm{nl}}$ | $\mathcal{N}_{x}\times\mathcal{N}_{y}$ | $\tau=0.002L$ | $\tau=0.003L$ | $l^{\textrm{nl}}$ | $\mathcal{N}_{x}\times\mathcal{N}_{y}$ | $\alpha=0.6$ | $\alpha=0.3$
$\overline{u}_{x}\times 100$ |  $\Delta$ [%] | $\overline{u}_{x}\times 100$ |  $\Delta$ [%] | $\overline{u}_{x}\times 100$ |  $\Delta$ [%] | $\overline{u}_{x}\times 100$ |  $\Delta$ [%]
$0.2L$ | $2\times 2$ | 6.12 | - | 6.13 | - | $0.25L$ | $2\times 2$ | 6.44 | - | 6.40 | -
$4\times 4$ | 6.27 | 2.59 | 6.27 | 2.32 | $4\times 4$ | 6.30 | 2.19 | 6.30 | 1.56
$5\times 5$ | 6.24 | 0.46 | 6.24 | 0.45 | $5\times 5$ | 6.29 | 0.15 | 6.29 | 0.20
$6\times 6$ | 6.25 | 0.24 | 6.25 | 0.25 | $6\times 6$ | 6.28 | 0.09 | 6.28 | 0.12
$0.4L$ | $2\times 2$ | 3.26 | - | 4.32 | - | $0.5L$ | $2\times 2$ | 6.76 | - | 6.92 | -
$4\times 4$ | 6.12 | 87.62 | 6.13 | 42.13 | $4\times 4$ | 6.39 | 5.54 | 6.42 | 7.27
$6\times 6$ | 6.22 | 1.79 | 6.23 | 1.72 | $5\times 5$ | 6.35 | 0.63 | 6.37 | 0.86
$8\times 8$ | 6.27 | 0.78 | 6.27 | 0.59 | $6\times 6$ | 6.32 | 0.35 | 6.34 | 0.47

#####  p-refinement

For a complete analysis of the convergence conditions derived from
h-refinement, we additionally perform a p-refinement fixing the parent mesh
parameters at $M^{\textrm{p}}=12$ and $r^{\textrm{p-c}}=1$. For this purpose,
$1^{\textrm{st}}$ and $2^{\textrm{nd}}$ order Lagrange shape functions are
chosen. Note that, here the p-refinement is only applicable to the parent
domain since the DOF (and hence the shape functions), following Eringen’s
approach are defined at the parent domain. The results of this analysis are
presented in Table 4. The nonlocal parameters assumed here are different from
that in h-refinement, to verify the universality of the convergence
conditions. As evident from Table 4, the relative change in the Euclidean
distance between the displacement fields obtained using the $1^{\textrm{st}}$
and $2^{\textrm{nd}}$ order Lagrange shape functions is less than the a priori
defined threshold (that is, $\Delta<2\%$), indicating converged simulations.

Table 4: Convergence of $\text{M}^{2}$-FEM solution to the nonlocal problem in the two case studies using p-refinement. Case study 1 | Case study 2
---|---
$l^{\textrm{nl}}$ | $\tau$ |  $1^{\textrm{st}}$ order |  $2^{\textrm{nd}}$ order | $\Delta$ | $l^{\textrm{nl}}$ | $\alpha$ |  $1^{\textrm{st}}$ order |  $2^{\textrm{nd}}$ order | $\Delta$
$\overline{u}_{x}\times 100$ | $\overline{u}_{x}\times 100$ | [%] | $\overline{u}_{x}\times 100$ | $\overline{u}_{x}\times 100$ | [%]
$0.2L$ | $0.002L$ | 6.19 | 6.15 | 0.65 | $0.3L$ | $0.6$ | 6.34 | 6.27 | 1.10
$0.003L$ | 6.26 | 6.28 | 0.32 | $0.3$ | 6.37 | 6.18 | 2.98
$0.3L$ | $0.002L$ | 6.08 | 6.15 | 1.15 | $0.5L$ | $0.6$ | 6.36 | 6.28 | 1.26
$0.003L$ | 6.21 | 6.29 | 1.29 | $0.3$ | 6.42 | 6.35 | 1.09

##  Impact of $\text{M}^{2}$-FEM on multiscale modeling

In order to highlight the significant impact of M2-FEM on multiscale modeling,
we showcase the ability of the underlying mesh-decoupling technique to: 1)
enable the use of unstructured meshes for the simulation of non-standard non-
rectangular geometries with (geometrically) non-similar nonlocal horizons; and
2) achieve high accuracy and high efficiency simultaneously, by tuning the
decoupled child and parent mesh-sizes, respectively.

###  Simulation of non-rectangular multiscale configurations

We consider a set of solids with non-rectangular shape and non-rectangular
nonlocal horizons. These choices enforce the use of unstructured mesh in both
the parent and child domains and hence, (recall from §1) cannot be simulated
by existing integral solvers. Three configurations, illustrated in Fig. (4)
along with their unstructured mesh, are considered. The three configurations
were chosen following a careful survey of some common geometrical features
underpinning multiscale solids in practical applications; they are not defined
arbitrarily to merely show the need for unstructured meshing. Further, the
nonlocal kernels used to simulate the configurations are also very closely
representative (in a functional sense) of the actual microstructural signature
in the emerging nonlocal effects noted in corresponding applications. For all
the configurations, the nonlocal horizon is truncated at points lying close to
and on the material boundaries [16, 15] resulting in non-constant irregular
shapes of the horizon. This phenomenon is schematically illustrated for
selected points in Fig. (4).

Figure 4: Schematic illustrating mesh and nonlocal kernel function
configurations in the three numerical test cases. a Nonlocal annulus: the
annular parent domain is discretized by unstructured triangular elements while
the child domain consisting of two orthogonal lines is discretized by 1D line
elements. b Nonlocal chip: the rectangular parent domain and the circular
child domain are discretized by structured quadrilateral elements and
unstructured triangular elements, respectively. c Nonlocal inclusion:
unstructured triangular mesh is used to discretize the whole domain and
adaptive mesh refinement technique is used in the nonlocal inclusion to
achieve higher accuracy. Different nonlocal kernel functions applied in the
three simulation cases are presented at the bottom-right corner of each
configuration.

####  Relevance and behavior of the selected configurations

The geometry, practical relevance, and mechanical behavior of the three
configurations are discussed below:

  * •

Nonlocal annulus: a nonlocal annular domain at the parent mesh level. At each
point within the annulus, the nonlocal horizon is assumed to be restricted
along bidirectional line segments. From a numerical perspective, the
differential dimensions assumed for the parent (2D) and child (1D) highlight
the significant flexibility afforded by $\text{M}^{2}$-FEM; indeed this choice
presents the most drastic case of geometric non-similarity between the parent
and child domains. Such reduced-order and preferential directions of
localization of nonlocal effects is typical of orthotropic and anisotropic
microstructures in composites and porous solids [5, 41]. The overall nonlocal
annulus geometry is commonly noted in naturally-occurring porous solids (such
as, animal bones and bamboo) and several real-world applications including
energy conversion devices (batteries) and biomedical implants [41, 2]. As
shown in [42, 8, 15], power-law kernels closely capture the deformation
behavior of these solids.

  * •

Nonlocal plane strain domain: a square-shaped parent domain whose material
points are associated with a circular nonlocal horizon (geometrical non-
similar to the square-shaped parent domain). This configuration could be very
relevant to nano- and micro- electromechanical systems [26], semiconductor
devices, and even periodic metastructures embedded with acoustic black holes
[43]. Exponential kernel functions are more representative of the dispersive
behavior in these structural configurations [26, 43].

  * •

Local plane strain domain with nonlocal inclusion: a square-shaped plane
strain local domain with an embedded circular nonlocal inclusion. The
inclusion is fully nonlocal, meaning that the nonlocal horizon at any material
point within the inclusion, coincides with the entire inclusion. This
configuration is an excellent example of interface problems and can be seen as
a multiscale extension of the classical inclusion problem and can be found
commonly in functionally-graded solids, layered media (such as, soil), and
solids with periodic reinforcements, interstitial defects, or distributed
damage [44, 32]. The strong nonlocality and dispersion noted in these
applications is well captured by rational kernels [44, 45].

##### Nonlocal annulus.

The inner and outer radii of the annulus are taken as
$r^{\textrm{in}}=0.3\textrm{m}$ and $r^{\textrm{out}}=0.5\textrm{m}$,
respectively. Further, the length scale of the nonlocal horizon is taken as
$l^{\textrm{nl}}=0.2\textrm{m}$. The parent and child domains are discretized
using triangular and linear Lagrange elements, respectively. The kernel
function is defined (similar to [8]) as:

$\displaystyle\mathcal{K}(x,y,x^{\prime},y^{\prime},\alpha)$
$\displaystyle=\frac{1}{2}\left[\mathcal{K}_{x}(x,y,x^{\prime},y^{\prime},\alpha)+\mathcal{K}_{y}(x,y,x^{\prime},y^{\prime},\alpha)\right]$
(53a) $\displaystyle\mathcal{K}_{x}(x,y,x^{\prime},y^{\prime},\alpha)$
$\displaystyle=\frac{1}{2\Gamma(1-\alpha)}\frac{1}{|x-x^{\prime}|^{\alpha}},\quad|x-x^{\prime}|<l^{\textrm{nl}},~{}y=y^{\prime}$
(53b) $\displaystyle\mathcal{K}_{y}(x,y,x^{\prime},y^{\prime},\alpha)$
$\displaystyle=\frac{1}{2\Gamma(1-\alpha)}\frac{1}{|y-y^{\prime}|^{\alpha}},\quad|y-y^{\prime}|<l^{\textrm{nl}},~{}x=x^{\prime}$
(53c)

where $\mathcal{K}_{x}$ and $\mathcal{K}_{y}$ are singular power-law functions
defined in the $x$ and $y$ directions. The above definition reduces the 2D
nonlocal horizon to bidirectional 1D line segments. The inner boundary of the
annulus is fixed, and the outer boundary is loaded with normal tractions of
constant magnitude $|\bm{t}^{\textrm{out}}|=0.1\textrm{Pa}$. The response of
the nonlocal annulus, in terms of the absolute value of the displacement
field, parameterized for different values of the fractional-order $\alpha$,
are presented in Fig. (5).

Figure 5: Response of the nonlocal annulus in terms of the absolute
displacement field, simulated via $\text{M}^{2}$-FEM. The results are
parameterized for different values of the nonlocal strength parameter
$\alpha$: a $\alpha=0.4$, b $\alpha=0.6$, c $\alpha=0.8$, and d $\alpha=1$.

As evident, $\text{M}^{2}$-FEM enables the smooth (lack of spurious gradients)
and non-divergent simulation of non-rectangular geometries. Observe that the
magnitude of the maximum absolute displacement (around the outer boundary)
increases as the power-law order $\alpha$ decreases. This simulated softening
behavior is consistent with the predictions of fractional-order approaches to
nonlocal elasticity [11]. Note that the power-law kernel in Eq. (53) is
denominated constant-order kernel in FC [30]. To this end, we highlight that
$\text{M}^{2}$-FEM can also be used to simulate variable- and distributed-
order kernels that lead to functionally-varying and functionally-distributed
space of the power-law exponent, respectively. These advanced kernels were
used recently to develop the highly generalized variable-order and
distributed-order nonlocal approaches [29, 11]. The asymmetric and anisotropic
kernels (unlike, Eringen’s approach [13]) used in these approaches enables the
simulation strongly nonlocal and multiscale effects induced by intricate
microstructures [20, 12].

##### Nonlocal plane strain domain.

This configuration, as shown in Fig. (4-b), consists of a square-shaped parent
domain with size $L=1\textrm{m}$, and a circular child domain with radius
$l^{\textrm{nl}}=0.2\textrm{m}$. The parent domain is discretized uniformly
using structured quadrilateral elements and the child domain is discretized
randomly using unstructured triangular elements. Note that the truncation of
nonlocal horizon near boundaries, requires the regeneration and re-
descretization of the corresponding child mesh. We adopt a 2D and radially
symmetric exponential kernel:

$\mathcal{K}(\bm{x},\bm{x^{\prime}})=\frac{1}{\tau_{1}}\mathrm{exp}\left(\frac{|\bm{x}-\bm{x}^{\prime}|}{\tau_{2}}\right),\quad|\bm{x}-\bm{x}^{\prime}|<l^{\textrm{nl}}$
(54)

where $\tau_{1}$ and $\tau_{2}$ are constitutive parameters that govern the
amplitude and decay rate of the kernel, respectively. Results of two
parametric studies within this example, for different values of $\tau_{1}$ and
$\tau_{2}$, are presented in Fig. (6) and Fig. (7), respectively. Aside from
the present discussion, an interesting result emerges from the two studies.
The results in Figs. (6,7) show that, depending on the value of parameters
$\tau_{1}$ and $\tau_{2}$, the nonlocal kernel in Eq. (54), can represent
either structural softening or stiffening behavior (a feature often difficult
to implement in multiscale constitutive modeling222A discussion on the
constitutive modeling of stiffening-softening effects is beyond the scope of
the present study. The interested reader is referred to [17, 37, 14, 46] for
detailed discussions on the ability of nonlocal constitutive theories in
modeling softening and/or stiffening effects.). More specifically, results in
Fig. (6-a) and Fig. (7-a) show the emergence of the softening behavior. The
softening effect is evident from the larger displacement amplitude (around the
central region of the domain) when compared with the local counterpart in Fig.
(6-d) and Fig. (7-d). On the other hand, a stiffening behavior results in
smaller displacement amplitudes (compared to the local response) as noted in
Figs. (6-c,7-c).

Figure 6: Response of the nonlocal plane strain domain simulated via
$\text{M}^{2}$-FEM, is presented in terms of the displacement fields along $x$
direction ($u_{x}$). The response is parameterized for different values of the
nonlocal strength $\tau_{2}$: a $\tau_{1}=1/1000$, b $\tau_{1}=1/2000$, and c
$\tau_{1}=1/3000$, with a fixed value of $\tau_{2}=1/100$. The local response
is presented in d for reference.

##### Local plane strain domain with nonlocal inclusion.

The outer dimensions of the local matrix and the radius of the inner inclusion
are assumed as $L_{x}\times L_{y}=1\textrm{m}\times 1\textrm{m}$ and
$r^{\textrm{in}}=0.15\textrm{m}$, respectively. Unstructured meshes with
triangular elements were used to discretize the matrix and the inclusion
(including both parent and child domains). Notably, mesh-decoupling directly
enabled the use of mesh-grading and mesh-refinement (see Fig. (8)) [39], to
reduce the overall computational cost. More specifically, these techniques
reduced the computational resources dedicated to the local matrix, while
simultaneously ensuring a detailed resolution of the nonlocal inclusion and
the interface. This result is of practical relevance in the simulation of more
complex configurations, such as those with distributed inclusions. Now, to
simulate this configuration we assume a rational kernel:

$\mathcal{K}(\bm{x},\bm{x^{\prime}})=\frac{1}{\tau_{1}|\bm{x}-\bm{x}^{\prime}|^{2}+\tau_{2}},\quad|\bm{x}-\bm{x}^{\prime}|\leq
2r^{\textrm{in}}$ (55)

where $\tau_{1}=100$ and $\tau_{2}=5000$ are nonlocal model parameters.
Displacement boundary conditions are employed:

$\bm{u}(x,y)\equiv\\{u_{x}(x,y),u_{y}(x,y)\\}=\\{x,0\\}\quad\forall(x,y)$ (56)

The simulated response is presented Fig. (8). Notice that the presence of
nonlocal effects distorts the displacement field within an appreciable
neighbourhood surrounding the interface. A clear change in displacement
gradient, also called ’weak discontinuity’, is observed along the
circumference of the interface. This behavior is typical of the response of
nonlocal inclusions [8, 44], and it is well captured by $\text{M}^{2}$-FEM.
The results are remarkable because, unlike existing approaches,
$\text{M}^{2}$-FEM readily discretizes the circular interface between the
nonlocal inclusion and the local matrix independently, while simultaneously
capturing the multiscale interaction via scale-bridging.

Figure 7: Response of the nonlocal plane strain domain simulated via
$\text{M}^{2}$-FEM, is presented in terms of the displacement fields along $x$
direction ($u_{x}$). The response is parameterized for different values of
nonlocal strength parameter $\tau_{2}$: a $\tau_{2}=1/100$, b $\tau_{2}=1/75$,
and c $\tau_{2}=1/50$, with a fixed value of $\tau_{1}=1/1000$. The local
response is presented in d for reference. Figure 8: Simulation of the local
plane strain domain with nonlocal inclusion and nonlocal interface, via
$\text{M}^{2}$-FEM. Results are presented in terms of the displacement field
along the: a $x$-direction ($u_{x}$), and b y-direction ($u_{y}$).
Orthographic views of the displacement fields are also presented to enable a
complete visualization.

###  High accuracy with high efficiency: overcoming restrictive trade-offs

The most significant impact of $\text{M}^{2}$-FEM, in view of the highly
multiscale character of emerging [20, 22] and future applications [7, 6], is
the ability to simultaneously achieve high accuracy and high efficiency while
requesting limited computational resource. This is in sharp contrast to
traditional solvers that impose implicit trade-offs between accuracy and
efficiency. Note that, differently from the approximation of local elastic
potential energies, multiscale constitutive models require the approximation
of additional integrals at each Gauss point in the parent domain, which
drastically increases time complexity of these nonlocal problems333Note that
the same concept also extends to gradient theories where additional gradients
would necessitate numerical evaluation. Although beyond the scope of this
study, recall from §1 that, gradient models are less general when compared to
integral models. In regards to the present discussion, they also lead to
practical difficulties in the numerical implementation, rooted in the
continuity of evaluated higher-order gradients.. As explained in detail in §1,
sharing the mesh between the parent and child domains drastically increases
the overall computational cost when tuning the element size mesh for better
accuracy in the nonlocal integrals444As seen in §4.2 and literature [27, 25],
the numerical accuracy of nonlocal finite element approaches is more sensitive
to the approximation of nonlocal integrals associated with the fine material
scales.. On the other hand, the decoupling of the meshes in
$\text{M}^{2}$-FEM, allows an independent tuning of the child element size for
better accuracy (via better approximation of the fine-scale nonlocal
integrals) and the parent element size for better efficiency. To enable a more
transparent visualization of this concept, we analyze the time complexity of
both $\text{M}^{2}$-FEM and traditional nonlocal finite element approaches
when simulating integral elasticity problems.

We consider a benchmark problem similar to that shown in Fig. (1). More
specifically, we consider a square parent domain (size $L\times L$) which is
discretized uniformly by the parent mesh with mesh number $M^{\textrm{p}}$ in
both the $x$ and $y$ directions. Further, the child domain at each parent
Gauss point is also assumed to be square-shaped (size $l^{\textrm{nl}}\times
l^{\textrm{nl}}$), and it is discretized uniformly by the child mesh with mesh
number $M^{\textrm{c}}$ in both the $x$ and $y$ directions. To simplify the
overall analysis, we do not consider the truncation of the child domain caused
by the presence of the material boundary and hence, ignore the change of child
mesh number for the parent Gauss points lying close to the boundaries. For the
above defined configuration, we obtain the total number of parent Gauss points
in the parent domain and the total number of child Gauss points in each child
domain, as:

$\displaystyle M^{g^{\textrm{p}}}$
$\displaystyle=g^{\textrm{p}}\left(M^{\textrm{p}}\right)^{2}$ (57a)
$\displaystyle M^{g^{\textrm{c}}}$
$\displaystyle=g^{\textrm{c}}\left(M^{\textrm{c}}\right)^{2}$ (57b)

where $g^{\textrm{p}}$ and $g^{\textrm{c}}$ are the number of parent and child
Gauss points in each parent and child element, respectively. Now, the total
number of individual evaluations and summations performed in the approximation
of the nonlocal integration and the weak-form integration, both following the
numerical quadrature rules, can be approximated (accurate to the limit of
multiplicative or additive constants) as:

$M^{\textrm{sum}}\approx
g^{\textrm{p}}g^{\textrm{c}}\left(M^{\textrm{c}}\right)^{2}\left(M^{\textrm{p}}\right)^{2}$
(58)

Note that the above expression is valid for all existing nonlocal finite
element approaches (including, $\text{M}^{2}$-FEM) for a first-order (one)
parent-child nesting; see similar analysis in [28]. With the above
understanding, recall the relative resolution $r^{\textrm{p-c}}$ of
$\text{M}^{2}$-FEM, used in §4.1 to represent the ratio between the parent and
child mesh size. By using the definitions of the relative resolution, and the
parent and child mesh numbers, $M^{\textrm{sum}}$ can be expressed as:

$M^{\textrm{sum}}\approx
g^{\textrm{p}}g^{\textrm{c}}\left(r^{\textrm{p-c}}\right)^{2}\left(\frac{l^{\textrm{nl}}}{L}\right)^{2}\left(M^{\textrm{p}}\right)^{4}\equiv
g^{\textrm{p}}g^{\textrm{c}}\frac{1}{\left(r^{\textrm{p-c}}\right)^{2}}\left(\frac{L}{l^{\textrm{nl}}}\right)^{2}\left(M^{\textrm{c}}\right)^{4}$
(59)

It immediately follows that, for a given configuration, the computational
complexity: 1) increases proportionally to the square of $r^{\textrm{p-c}}$
when fixing the parent mesh number $M^{\textrm{p}}$, and 2) decreases with the
square of $r^{\textrm{p-c}}$ when fixing the child mesh number
$M^{\textrm{c}}$. These remarkable characteristics are valid only for
$\text{M}^{2}$-FEM and are not applicable to traditional nonlocal finite
element approaches because, the relative resolution is fixed as
$r^{\textrm{p-c}}=1$ following the direct sharing of the mesh between the
parent and child domains. To complement the above theoretical deduction, we
present performance data from simulations in the Fig. (9).

Figure 9: Variation of both computational time and error with respect to the
parent and child mesh numbers. a 3D bar plot showing growth in computational
time following an increase in both the parent and child mesh numbers. b
Variation of computational time cost (left $y$-axis and blue curves) and the
relative error (right $y$-axis and orange curves) with different mesh numbers.
The legend specifying different line styles applies to both blue and orange
curves.

Figure (9) compares the variation of the computational cost (in terms of the
average time needed for a set of simulations) against the relative error (in
terms of the difference in simulated $\textrm{max}(u_{x})$ from the reference
value) for Case Study 1. All the simulations were performed on a personal
computer equipped with an Intel(R) Core(TM) i7-9750H processor and 16 GB RAM.
Figure (9-a) presents the growth of total computational time with respect to
both $M^{\textrm{p}}$ and $r^{\textrm{p-c}}$. As evident from the expression
of $M^{\textrm{sum}}$ in Eq. (59), the time complexity can be controlled by
varying either the parent ($M^{\textrm{p}}$) or child mesh size
($r^{\textrm{p-c}}$). Figure (9-b) presents the variation of time cost and
error with respect to $M^{\textrm{p}}$ and $r^{\textrm{p-c}}$. As evident, by
tuning the values of $M^{\textrm{p}}$ and $r^{\textrm{p-c}}$,
$\text{M}^{2}$-FEM can be highly efficient (reflect in its ability to save
computational time) while still maintaining the same level of accuracy. This
aspect is more transparent from two sets of simulations in the Fig. (9):

  * •

the solution obtained at $(M^{\textrm{p}},r^{\textrm{p-c}})=(10,2)$ requires
(roughly) the same simulation time as that at
$(M^{\textrm{p}},r^{\textrm{p-c}})=(15,1)$, but has much smaller error (see
the data highlighted by dotted cyan circles),

  * •

the solution obtained at $(M^{\textrm{p}},r^{\textrm{p-c}})=(10,3)$ has
roughly the same error as that at $(M^{\textrm{p}},r^{\textrm{p-c}})=(25,1)$,
but requires less computational time (see the data highlighted by dotted
magenta circles).

Notably, for a fixed parent mesh, the simulation using a finer child mesh
($r^{\textrm{p-c}}>1$) shows better performance than the one with a non-
refined child mesh ($r^{\textrm{p-c}}=1$). These results indicate that, when
compared to existing (mesh-coupling-based) integral solvers,
$\text{M}^{2}$-FEM equipped with mesh-decoupling technique admits greater
computational flexibility; ultimately, leading to higher accuracy and higher
efficiency within fixed computational resources.

##  Summary and outlook

This study presented the formulation, the numerical implementation, and the
performance of $\text{M}^{2}$-FEM that is an advanced finite element solver
capable of simulating generalized integral theories. $\text{M}^{2}$-FEM
derives its unique characteristics from the concept of multiscale mesh-
decoupling and scale-bridging, which combined create a numerically independent
and localized set of meshes to approximate the underlying integral theory at
different material scales, in a physically consistent manner. The
characteristic multiscale independence has a multifold impact on the overall
computational solver. More specifically, it enables the independent topology-
tuning of the scale-specific localized meshes such that $\text{M}^{2}$-FEM is
able to simulate complex multiscale configurations and admit generalized
kernel functions. These latter characteristics are achieved while
simultaneously providing high accuracy and computational efficiency. These
features enable $\text{M}^{2}$-FEM to model the static and dynamic behavior of
a wide class of advanced materials that exhibit unconventional geometry,
nonlocal and multiscale effects, and anomalous dispersion behavior.

Overall, $\text{M}^{2}$-FEM removes significant computational drawbacks in the
application of integral structural theories. While the present study was
developed based on Eringen’s strain-driven nonlocal theory, $\text{M}^{2}$-FEM
is very general in nature and can be readily applied to other integral
theories including, for example, self-consistent approaches, homogenization
approaches, multiscale averaging approaches, and displacement-driven nonlocal
approaches. The extension of $\text{M}^{2}$-FEM to these approaches merely
requires a different numerical evaluation of the child-domain integral
contributions by following the model-specific constitutive theory. The
underlying mesh-decoupling and scale-bridging, fundamental pillars of
$\text{M}^{2}$-FEM, remain unaltered. The recursive nesting ability of
$\text{M}^{2}$-FEM highlighted the ability of the proposed method to support
the analysis of advanced material models with hierarchical nesting. These
characteristics suggest that $\text{M}^{2}$-FEM can serve as the foundation
for a next-generation class of solvers offering critical capabilities to
enable the simulation of highly multiscale systems.

##### Data Availability.

All the necessary data and information required to reproduce the results are
available in the paper and the supplementary information document.

##### Declaration of Competing Interest.

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##### Acknowledgements.

The authors gratefully acknowledge the financial support of the National
Science Foundation under grants MOMS #1761423 and DCSD #1825837. S.P.
acknowledges the partial support received through Dimitris N. Chorafas
Foundation Award (2022). Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the author(s) and do
not necessarily reflect the views of the National Science Foundation.

## References

  * [1] Ganti, S. S., Liu, T.-W. & Semperlotti, F.  Topological edge states in phononic plates with embedded acoustic black holes.  _Journal of Sound and Vibration_ 466, 115060 (2020).
  * [2] Mannan, S., Paul Knox, J. & Basu, S.  Correlations between axial stiffness and microstructure of a species of bamboo.  _Royal Society Open Science_ 4, 160412 (2017).
  * [3] Parameswaran, V. & Shukla, A.  Processing and characterization of a model functionally gradient material.  _Journal of Materials Science_ 35, 21–29 (2000).
  * [4] Russillo, A. F. & Failla, G.  Wave propagation in stress-driven nonlocal Rayleigh beam lattices.  _International Journal of Mechanical Sciences_ 215, 106901 (2022).
  * [5] Kanouté, P., Boso, D., Chaboche, J.-L. & Schrefler, B.  Multiscale methods for composites: a review.  _Archives of Computational Methods in Engineering_ 16, 31–75 (2009).
  * [6] Hoekstra, A., Chopard, B. & Coveney, P.  Multiscale modelling and simulation: a position paper.  _Philosophical Transactions of the Royal Society A_ 372, 20130377 (2014).
  * [7] Van Der Giessen, E. _et al._ Roadmap on multiscale materials modeling.  _Modelling and Simulation in Materials Science and Engineering_ 28, 043001 (2020).
  * [8] Patnaik, S., Jokar, M., Ding, W. & Semperlotti, F.  On the role of the microstructure in the deformation of porous solids.  _npj Computational Materials_ 8, 1–7 (2022).
  * [9] Mei, C. _et al._ A nonlocality-based homogenization method for dynamics of metamaterials.  _Composite Structures_ 295, 115716 (2022).
  * [10] Madeo, A., Barbagallo, G., d’Agostino, M. V., Placidi, L. & Neff, P.  First evidence of non-locality in real band-gap metamaterials: determining parameters in the relaxed micromorphic model.  _Proceedings of the Royal Society A_ 472, 20160169 (2016).
  * [11] Ding, W., Patnaik, S. & Semperlotti, F.  Multiscale nonlocal elasticity: A distributed order fractional formulation.  _International Journal of Mechanical Sciences_ 226, 107381 (2022).
  * [12] Ding, W., Patnaik, S. & Semperlotti, F.  Multiscale nonlocal beam theory: An application of distributed-order fractional operators (2022).  Preprint available at https://arxiv.org/abs/2207.13787.
  * [13] Batra, R.  Misuse of Eringen’s nonlocal elasticity theory for functionally graded materials.  _International Journal of Engineering Science_ 159, 103425 (2021).
  * [14] Patnaik, S., Sidhardh, S. & Semperlotti, F.  Displacement-driven approach to nonlocal elasticity.  _European Journal of Mechanics-A/Solids_ 92, 104434 (2022).
  * [15] Patnaik, S. & Semperlotti, F.  A generalized fractional-order elastodynamic theory for non-local attenuating media.  _Proceedings of the Royal Society A_ 476, 20200200 (2020).
  * [16] Polizzotto, C.  Nonlocal elasticity and related variational principles.  _International Journal of Solids and Structures_ 38, 7359–7380 (2001).
  * [17] Eringen, A. C.  Theory of nonlocal elasticity and some applications.  Tech. Rep., Princeton Univ NJ Dept of Civil Engineering (1984).
  * [18] McDowell, D. L.  A perspective on trends in multiscale plasticity.  _International Journal of Plasticity_ 26, 1280–1309 (2010).
  * [19] Liu, S., Duan, K., Li, L., Wang, X. & Hu, Y.  A multilayer coarse-grained molecular dynamics model for mechanical analysis of mesoscale graphene structures.  _Carbon_ 178, 528–539 (2021).
  * [20] Patnaik, S., Jokar, M., Ding, W. & Semperlotti, F.  Distillation of nonlocality in porous solids (2022).  Preprint available at https://www.researchsquare.com/article/rs-1875929/v1.
  * [21] Stempin, P. & Sumelka, W.  Formulation and experimental validation of space-fractional Timoshenko beam model with functionally graded materials effects.  _Computational Mechanics_ 68, 697–708 (2021).
  * [22] Failla, G. & Zingales, M.  Advanced materials modelling via fractional calculus: challenges and perspectives.  _Philosophical Transactions of the Royal Society A_ 378, 20200050 (2020).
  * [23] Nair, S.  _Nonlocal Acoustic Black Hole Metastructures: Achieving Ultralow Frequency and Broadband Vibration Attenuation_.  Ph.D. thesis, Purdue University Graduate School (2019).
  * [24] Chopard, B., Borgdorff, J. & Hoekstra, A. G.  A framework for multi-scale modelling.  _Philosophical Transactions of the Royal Society A_ 372, 20130378 (2014).
  * [25] Pisano, A., Sofi, A. & Fuschi, P.  Nonlocal integral elasticity: 2D finite element based solutions.  _International Journal of Solids and Structures_ 46, 3836–3849 (2009).
  * [26] Sidhardh, S. & Ray, M.  Effect of nonlocal elasticity on the performance of a flexoelectric layer as a distributed actuator of nanobeams.  _International Journal of Mechanics and Materials in Design_ 14, 297–311 (2018).
  * [27] Patnaik, S., Sidhardh, S. & Semperlotti, F.  A Ritz-based finite element method for a fractional-order boundary value problem of nonlocal elasticity.  _International Journal of Solids and Structures_ 202, 398–417 (2020).
  * [28] Patnaik, S., Sidhardh, S. & Semperlotti, F.  Fractional-order models for the static and dynamic analysis of nonlocal plates.  _Communications in Nonlinear Science and Numerical Simulation_ 95, 105601 (2021).
  * [29] Patnaik, S., Jokar, M. & Semperlotti, F.  Variable-order approach to nonlocal elasticity: Theoretical formulation, order identification via deep learning, and applications.  _Computational Mechanics_ 69, 267–298 (2021).
  * [30] Das, S.  _Functional fractional calculus_ , vol. 1 (Springer, 2011).
  * [31] Ostoja-Starzewski, M.  Microstructural randomness versus representative volume element in thermomechanics.  _Journal of Applied Mechanics_ 69, 25–35 (2002).
  * [32] Romkes, A. & Oden, J. T.  Adaptive modeling of wave propagation in heterogeneous elastic solids.  _Computer Methods in Applied Mechanics and Engineering_ 193, 539–559 (2004).
  * [33] You, H., Yu, Y., Trask, N., Gulian, M. & D’Elia, M.  Data-driven learning of nonlocal physics from high-fidelity synthetic data.  _Computer Methods in Applied Mechanics and Engineering_ 374, 113553 (2021).
  * [34] Flaschel, M., Kumar, S. & De Lorenzis, L.  Unsupervised discovery of interpretable hyperelastic constitutive laws.  _Computer Methods in Applied Mechanics and Engineering_ 381, 113852 (2021).
  * [35] Gallet, A. _et al._ Structural engineering from an inverse problems perspective.  _Proceedings of the Royal Society A_ 478, 20210526 (2022).
  * [36] Jokar, M. & Semperlotti, F.  Finite element network analysis: A machine learning based computational framework for the simulation of physical systems.  _Computers & Structures_ 247, 106484 (2021).