\section{Introduction}

The failure mechanism of structure components and solid materials is fundamental for the research of structural integrity. Classical continuum mechanics (CCM) employs spatial derivatives to describe the material behaviors with a requirement of continuum material properties during the deformation. However, derivatives are undefined at discontinuities of materials. Thus the CCM is inherently difficult to material failure. In 2000, the peridynamics (PD) theory \citep{silling2000reformulation} was introduced to remove this drawback. The critical idea of the PD theory is that the PD employs spatial integrals rather than derivatives. Since the integrals are still defined at discontinuities, the PD is suitable for describing non-continuum materials.

The first proposed PD theory is the bond-based formulation \citep{silling2000reformulation} which treats the bond between two material points as a spring, resulting in a restriction of a fixed Poisson's ratio of 1/4 \citep{gerstle2005peridynamic,madenci_peridynamic_2014}. This substantial limitation blocks the PD from a broader application with various materials. Thus, the fixed Poisson's ratio requirement must be removed. Indeed, \citet{silling_peridynamic_2007} extended the PD and developed the state-based PD (SBPD) to eliminate this constraint. In the SBPD, the bond force is redefined not only by the bond between the two material points but also by other neighboring bonds. \citet{silling_peridynamic_2007} provided two versions of the SBPD: the ordinary state-based PD (OSBPD) and the non-ordinary state-based PD (NOSBPD). The OSBPD divided the stretch of a bond into a hydrostatic and a deviatoric part, similar to dividing the strain tensor into hydrostatic and deviatoric strain tensors. The NOSBPD introduced a non-local displacement gradient and employed classical constitutive law to describe the bond forces. Both OSBPD and NOSBPD can remove the restriction of Poisson's ratio. Note that NOSBPD is inherent unstable, thus requires zero-energy mode control \cite{breitenfeld_non-ordinary_2014}. Thus, in this paper, an extended OSBPD is developed based on the OSBPD model.

Although the SBPD provides a theoretical framework for the discontinuity analysis with various materials, there are still three major issues that attract researchers to improve the generosity and decrease the computational cost: (a) the surface correction, (b) the volume correction, and (c) the meshless method with uniform discretization. All three problems are related to the spherical horizons in the PD. The PD horizon is defined as an interaction domain where the classical strain energy equals that of a PD. The requirement of a spherical horizon might be reasonable in the inner areas away from the structure boundaries. Nevertheless, as we get close to the edges of the structure, the spherical horizons are usually truncated, called the surface effect in PD. Several methods have been proposed to reduce or eliminate the surface effect. For example, \citet{gerstle2005peridynamic} and \citet{oterkus_peridynamic_2014} introduced fictitious nodes outsides of the boundary which may only be applied for simple geometries. \citet{scabbia2021novel} proposed a revised fictitious node approach with an extrapolation by Taylor series expansion to reduce the surface effect. \citet{madenci_peridynamic_2014} computed approximate corrections of the material constants for nodes near the surface by equating local SED to the non-local SED. This approach is exact in only homogeneous deformations because it assumes six simple loading conditions applied on the whole domain for surface correction:
\begin{subequations}\label{eq:loadcondi}
\begin{align}
    \varepsilon_{11}\neq 0,\ \varepsilon_{22}=\varepsilon_{33}=\varepsilon_{12}=\varepsilon_{23}=\varepsilon_{31}=0,\\
     \varepsilon_{22}\neq 0,\ \varepsilon_{11}=\varepsilon_{33}=\varepsilon_{12}=\varepsilon_{23}=\varepsilon_{31}=0,\\
      \varepsilon_{33}\neq 0,\ \varepsilon_{11}=\varepsilon_{22}=\varepsilon_{12}=\varepsilon_{23}=\varepsilon_{31}=0,\\
       \varepsilon_{12}\neq 0,\ \varepsilon_{11}=\varepsilon_{22}=\varepsilon_{33}=\varepsilon_{23}=\varepsilon_{31}=0,\\
        \varepsilon_{23}\neq 0,\ \varepsilon_{11}=\varepsilon_{22}=\varepsilon_{33}=\varepsilon_{12}=\varepsilon_{31}=0,\\
         \varepsilon_{31}\neq 0,\ \varepsilon_{11}=\varepsilon_{22}=\varepsilon_{33}=\varepsilon_{12}=\varepsilon_{23}=0.
\end{align}
\end{subequations}
\citet{shen_peridynamic_2021} proposed a similar surface correction method yet requires the FEM simulation performed forward in detail.

Moreover, the numerical implementation of the PD is generally by a simple meshfree method with a one-point integration method \citep{silling_meshfree_2005}. Such quadrature requires a volume correction to improve the integration accuracy. The simplest volume correction method is introducing a linear correction factor \citep{silling_meshfree_2005,parks_implementing_2008,le2014two}. Some other more accurate correction methods can be found in \citep{yu_new_2011,seleson_improved_2014,ren_dual-horizon_2017,seleson_numerical_2018}, but usually with a higher computational cost.

Furthermore, because of the spherical horizon restriction, one has to use regular uniform numerical discretization with constant horizon size to solve the PD equations instead of irregular non-uniform meshes with various horizon sizes. Otherwise, one may run into a "ghost force" problem which breaks the conservation of linear and angular momentum. However, it is also clear that uniform discretization has two drawbacks: (a). A high computational cost because the refinement level is applied anywhere in the domain, (b). The uniform mesh is difficult to conformal to complex geometries. To solve the "ghost force" issue and to allow numerical implementation of non-uniform discretization with various horizon sizes, \citet{ren_dual-horizon_2016,ren_dual-horizon_2017} proposed a dual-horizon PD model.

To solve the three critical problems of volume and surface corrections and uniform discretization implementation induced by the spherical horizon as a whole, many researchers have proposed a few PD models with non-spherical horizons. For example, \citet{madenci_peridynamic_2019} developed a PD-least square minimization (PDLSM), and \citet{liu_revised_2021} formulated a revised NOSBPD based on Taylor series expansion and least square minimization, which works for arbitrary horizon shapes. Still, there are instability problems in the dynamical analysis of the PDLSM and NOSBPD. \citet{madenci_bond-based_2021} derived bond-based peridynamics (BBPD) with stretch and rotation Kinematics based on the PDLSM. \citet{hu_thermomechanical_2018} presented a generalized BBPD and OSBPD models with non-spherical horizons, where they also assumed a simple loading condition applied to the whole domain to approximately correct material constants, i.e. $\varepsilon_{11}=\varepsilon_{22}=\varepsilon_{33}\neq 0$ and $\varepsilon_{12}=\varepsilon_{23}=\varepsilon_{31}=0$. \citet{mitchell_position-aware_2015} demonstrated a position-aware OSBPD model for non-spherical horizons in which the six simple loading conditions of \cref{eq:loadcondi} were assumed for correcting the scalar force states. Here, we must point out that the simple homogeneous deformation assumptions may bring errors because the real domain solutions are generally complex. Thus, an SBPD model is still needed without the assumptions given by \cref{eq:loadcondi}.

Based on the OSBPD, this paper presents an extended OSBPD (XOSBPD) model for non-spherical horizons by introducing the Lagrange multipliers to guarantee the non-local dilatation and non-local strain energy density (SED) are equal to local dilatation and local SED, respectively. The significant advantages of the XOSBPD are: (a). the XOSBPD does not assume the simple loading condition given in \cref{eq:loadcondi} and works for arbitrary deformation conditions, (b). no surface or volume corrections are needed, and (c). the XOSBPD allows non-uniform discretization implementation with various horizon sizes.

This paper is organized as follows. The OSBPD theory is briefly reviewed in \cref{sec:OSBPD}, and the XOSBPD is proposed and presented in \cref{sec:XOSBPD}. After that, the numerical implementation of XOSBPD is described in \cref{sec:NumImp}. Then, three classic examples are performed to demonstrate the proposed XOSBPD in \cref{sec:NumR}. Finally, the conclusions are drawn in \cref{sec:conc}.

\section{Ordinary state based PD}\label{sec:OSBPD}

In this section, the OSBPD for the 3-D domain first developed by \citet{silling_peridynamic_2007} is briefly reviewed for completeness. The equation of motion in the PD is given by:
\begin{equation}\label{eq:eqm}
    \rho \Ddot{\textbf{u}}=\textbf{L}(\textbf{x},t)+\textbf{b}(\textbf{x},t),
\end{equation}
where $\rho$ is the material mass density, $\textbf{u}$ is the displacement vector field, $\textbf{L}$ is non-local internal force density, and $\textbf{b}$ is externally applied body force density. The $\textbf{L}$ can be evaluated as:
\begin{equation}\label{eq:sillL}
   \textbf{L}(\textbf{x},t)=\int_{H_x} \left\{ \underline{\textbf{T}}(\textbf{x}) \langle \boldsymbol{\xi} \rangle - \underline{\textbf{T}}(\textbf{x}') \langle -\boldsymbol{\xi} \rangle \right\} \mathrm{d} V_{x'},
\end{equation}
where $H_x$ is the horizon of point $\textbf{x}$ and is required to be a sphere for both BBPD and SBPD \cite{silling2000reformulation,silling_peridynamic_2007}. $\boldsymbol{\xi}$ is a bond vector defined as:
\begin{equation}
    \boldsymbol{\xi}=\textbf{x}'-\textbf{x},
\end{equation}
and $\underline{\textbf{T}}(\textbf{x}) \langle \boldsymbol{\xi} \rangle$ is the force state. The \emph{state} is a mathematical object represented with a underline which maps the bond vector $\boldsymbol{\xi}$ to a tensor. For example, the un-deformed bond length state $ \underline{x}$ is defined as:
\begin{equation}
    \underline{x}\langle \boldsymbol{\xi} \rangle=|\boldsymbol{\xi}|.
\end{equation}
And the deformation vector state $\underline{\boldsymbol{Y}}$ is defined as:
\begin{equation}
    \underline{\boldsymbol{Y}}\langle \boldsymbol{\xi} \rangle=\boldsymbol{\xi}+\boldsymbol{\eta},
\end{equation}
in which $\boldsymbol{\eta}=\textbf{u}\left(\textbf{x}'\right)-\textbf{u}\left(\textbf{x}\right)$ is the relative displacement of the bond. One important scalar state in OSBPD is the extension state $\underline{e}$ defined as:
\begin{equation}
    \underline{e}\langle \boldsymbol{\xi} \rangle=|\underline{\boldsymbol{Y}}|-\underline{x}.
\end{equation}
Similar to dividing the strain into the hydrostatic and deviatoric parts in the CCM, the extension state can also divided into two parts:
\begin{equation}
    \underline{e}=\underline{e}^h+\underline{e}^d,
\end{equation}
in which the superscripts $h$ and $d$ represent the hydrostatic and deviatoric, respectively. The hydrostatic extension state is defined as:
\begin{equation}
    \underline{e}^h= \frac{\theta \underline{x}}{ 3},
\end{equation}
where $\theta$ is the non-local dilatation and for a spherical horizon:
\begin{equation}\label{eq:theta_sph}
    \theta =\theta \left(\underline{e} \right) =\frac{3}{m}\left(\underline{\omega} \underline{x} \right)\bullet \underline{e},
\end{equation}
where $\underline{\omega}\langle \boldsymbol{\xi} \rangle=\omega\left(|\boldsymbol{\xi}| \right)$ is an influence function and $m$ is the weighted volume, defined as:
\begin{equation}
  m=  \left(\underline{\omega} \underline{x} \right)\bullet \underline{x}.
\end{equation}
Here, the dot product of two scalar states of $\underline{a}$ and $\underline{b}$ is defined as:
\begin{equation}
    \underline{a}  \bullet \underline{b}=\int_{H_x} \underline{a}\langle \boldsymbol{\xi} \rangle \underline{b}\langle \boldsymbol{\xi} \rangle \mathrm{d} V_{x'}.
\end{equation}
For example, the weighted volume is evaluated as:
\begin{equation}
    m=\int_{H_x} \omega\left(|\boldsymbol{\xi}| \right) |\boldsymbol{\xi}|^2 \mathrm{d} V_{x'}=\int_{0}^{\delta} \int^{2\pi}_0 \int^{\pi}_0  \omega(r) r^4 \sin \phi \mathrm{d} \phi \mathrm{d} \theta \mathrm{d}r=4\pi \int_{0}^{\delta} \omega(r)r^4 \mathrm{d}r,
\end{equation}
where $\delta$ is the horizon size.

For the OSBPD, the force vector states $\underline{\textbf{T}}(\textbf{x}) \langle \boldsymbol{\xi} \rangle$ is defined as:
\begin{equation}\label{eq:T}
    \underline{\textbf{T}}(\textbf{x}) \langle \boldsymbol{\xi} \rangle=\underline{t}(\textbf{x}) \langle \boldsymbol{\xi} \rangle \frac{\boldsymbol{\eta}+\boldsymbol{\xi}}{|\boldsymbol{\eta}+\boldsymbol{\xi}|},
\end{equation}
where $\underline{t}(\textbf{x}) \langle \boldsymbol{\xi} \rangle$ is the scalar force state and is derived from the Fréchet derivative of the strain energy density with respect to the extension state as:
\begin{equation}\label{eq:frederi}
    \underline{t}(\textbf{x}) \langle \boldsymbol{\xi} \rangle =\left(\frac{\partial W_{pd}}{\partial \underline{e}} \right)^f =\left(\frac{\partial W_{pd}}{\partial \underline{e}^h} \right)^f + \left(\frac{\partial W_{pd}}{\partial \underline{e}^d} \right)^f,
\end{equation}
in which the superscript $f$ represents the Fréchet derivative is defined in \citep{silling_peridynamic_2007}.
For an isotropic material and a spherical horizon, the non-local strain energy density can be expressed as:
\begin{equation}\label{eq:sed_sph}
    W_{pd}=\frac{\kappa \theta^2}{2}+\frac{15}{2m}\mu(\underline{\omega} \underline{e}^d )\bullet \underline{e}^d,
\end{equation}
in which $\kappa$ is the bulk modulus and $\mu$ is the shear modulus. Substituting \cref{eq:sed_sph} into \cref{eq:frederi}, we have the scalar force state:
\begin{equation}\label{eq:finscalfor_sph}
     \underline{t}(\textbf{x}) \langle \boldsymbol{\xi} \rangle=\frac{3 \kappa \theta}{m} \underline{\omega} \underline{x} + \frac{15}{m} \mu \underline{\omega} \underline{e}^d.
\end{equation}
It is worth noting that \cref{eq:theta_sph,eq:sed_sph} are only valid for the spherical horizon. Thus, \cref{eq:finscalfor_sph} only works for the spherical horizon.

\section{Extended ordinary state based PD}\label{sec:XOSBPD}

In this section, we extend the horizons of the OSBPD from sphere to non-sphere ones. Thus removing the requirements of volume and surface corrections and enabling non-uniform discretization implementation with various horizon sizes. The presented XOSBPD introduces two correction items for the non-spherical horizon and arbitrary deformation conditions.

As described in \cref{sec:OSBPD}, the essential idea of the OSBPD are: (a) dividing the extension state $\underline{e}$ into hydrostatic extension state $\underline{e}^h$ and deviatoric extension state $\underline{e}^d$, (b) expressing the non-local strain energy density in terms of the $\underline{e}^h$ and $\underline{e}^d$, and then (c) using \cref{eq:frederi} to derive the scalar force state $\underline{t}$. The splitting of the extension state $\underline{e}$ is based on the assumption of small deformation and the classical kinematics as:
\begin{equation}
\begin{split}
     \underline{e}&=|\underline{\boldsymbol{Y}}|-\underline{x}=\frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}\boldsymbol{\xi}}{|\boldsymbol{\xi}|}\\
     &=\frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}^h\boldsymbol{\xi}}{|\boldsymbol{\xi}|} + \frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}^d\boldsymbol{\xi}}{|\boldsymbol{\xi}|}\\
     &=\underline{e}^h+\underline{e}^d,
\end{split}
\end{equation}
in which $\boldsymbol{\varepsilon}$, $\boldsymbol{\varepsilon}^h$, and $\boldsymbol{\varepsilon}^d$ are the strain tensor, hydrostatic strain tensor, and deviatoric strain tensor, respectively. Note that the non-local dilatation must equal to the local dilatation, that is $\theta=\varepsilon_{ii}$. Thus, we have:
\begin{align}\label{eq:eh_xpd}
    \underline{e}^h & =\frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}^h\boldsymbol{\xi}}{|\boldsymbol{\xi}|}=\frac{\boldsymbol{\xi}\cdot \frac{\theta \textbf{I}}{n_d}\boldsymbol{\xi}}{|\boldsymbol{\xi}|}=\frac{\theta \underline{x}}{n_d},\\
    \label{eq:ed_xpd}
    \underline{e}^d & =\frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}^d\boldsymbol{\xi}}{|\boldsymbol{\xi}|} =\underline{e}- \frac{\theta \underline{x}}{n_d},
\end{align}
in which $n_d$ is the number of dimensions.
From the CCM, the strain energy density of an isotropic material is:
\begin{equation}\label{eq:sde_ccm}
    W_{CCM}= \frac{1}{2} \kappa \varepsilon_{ii}\varepsilon_{jj} + \mu \boldsymbol{\varepsilon}^d_{ij}\boldsymbol{\varepsilon}^d_{ij}.
\end{equation}
As stated in \cite{silling_peridynamic_2007}, the non-local SED defined in \cref{eq:sed_sph} equals to that of the local SED only for the spherical horizon. To extend to the OSBPD for non-spherical horizons, we redefined the non-local strain energy density as:
\begin{equation}\label{eq:sde_xpd}
     W_{pd}=\frac{\kappa \theta^2}{2}+\mu(\underline{\omega}^d \underline{e}^d )\bullet \underline{e}^d.
\end{equation}
In this study, we refer the $\underline{\omega}^d\langle \boldsymbol{\xi} \rangle=\omega^d( \boldsymbol{\xi})$ as deviatoric influence function. For spherical horizon, ${\omega}^d$ is defined as:
\begin{equation}
   {\omega}^d={\omega}^d_s,
\end{equation}
in which ${\omega}^d_s$ is the spherical deviatoric influence function, and comparing \cref{eq:sed_sph} with \cref{eq:sde_xpd} reveals that the ${\omega}^d_s$ is:
\begin{equation}
    {\omega}^d_s=\frac{n_d(n_d+2)}{2m}\omega\left(|\boldsymbol{\xi}| \right),
\end{equation}
where the weighted volume $m$ is defined as:
\begin{equation}
    m=
    \begin{cases}
     &4\pi \int_{0}^{\delta} \omega(r)r^4 \mathrm{d}r, \quad \text{3-D},\\
     &2\pi \int_{0}^{\delta} \omega(r)r^3 \mathrm{d}r, \quad \text{2-D}.
    \end{cases}
\end{equation}
Besides, the definition of the non-local dilatation in \cref{eq:theta_sph} can ensure $\theta=\varepsilon_{ii}$ only for spherical horizon. For non-spherical horizon, we redefine the non-local dilatation as below:
\begin{equation}\label{eq:theta_xpd}
    \theta=(\underline{\omega}^h \underline{x}) \bullet \underline{e}.
\end{equation}
In this work, we refer the $\underline{\omega}^h\langle \boldsymbol{\xi} \rangle=\omega^h( \boldsymbol{\xi})$ as hydrostatic influence function. For spherical horizon, ${\omega}^h$ is defined as:
\begin{equation}
     {\omega}^h={\omega}^h_s,
\end{equation}
in which ${\omega}^h_s$ is the spherical hydrostatic influence function. Comparing \cref{eq:theta_sph} and \cref{eq:theta_xpd} reveals that the ${\omega}^h_s$ is expressed as:
\begin{equation}
    {\omega}^h_s=\frac{n_d}{m}\omega\left(|\boldsymbol{\xi}| \right).
\end{equation}
Equating \cref{eq:sde_xpd} to \cref{eq:sde_ccm} leads to:
\begin{equation}\label{eq:condi}
    \theta=\varepsilon_{ii}\quad  \mathrm{and}\quad (\underline{\omega}^d \underline{e}^d )\bullet \underline{e}^d=\boldsymbol{\varepsilon}^d:\boldsymbol{\varepsilon}^d.
\end{equation}
In the following section, we derive the hydrostatic and the deviatoric influence functions, satisfying the \cref{eq:condi} for arbitrary horizon shapes.

\subsection{Hydrostatic influence function}

Consider a 3-D body ($n_d=3$), we rearrange $\theta=\varepsilon_{ii}$ as:
\begin{equation}\label{eq:derivOmegaH}
\begin{split}
   0&=(\underline{\omega}^h \underline{x}) \bullet \underline{e}- \left(\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33}\right)\\
   &=\int_{H_x} {\omega}^h |\boldsymbol{\xi}| e \mathrm{d}V_{x'} - \left(\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33}\right)\\
    &=\int_{H_x} {\omega}^h |\boldsymbol{\xi}| \frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}\boldsymbol{\xi}}{|\boldsymbol{\xi}|} \mathrm{d}V_{x'}- \left(\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33}\right)\\
    &=\varepsilon_{11}\left( \int_{H_x} {\omega}^h \xi_1^2 \mathrm{d}V_{x'} -1 \right) + \varepsilon_{22}\left(\int_{H_x} {\omega}^h \xi_2^2\mathrm{d}V_{x'} -1 \right) + \varepsilon_{33}\left(\int_{H_x} {\omega}^h \xi_3^2\mathrm{d}V_{x'} -1 \right)\\
    &+\varepsilon_{12}\int_{H_x} {\omega}^h \xi_1\xi_2\mathrm{d}V_{x'} + \varepsilon_{23}\int_{H_x} {\omega}^h \xi_2\xi_3\mathrm{d}V_{x'} + \varepsilon_{31}\int_{H_x} {\omega}^h \xi_3\xi_1\mathrm{d}V_{x'}.
\end{split}
\end{equation}
For an arbitrary strain tensor $\boldsymbol{\varepsilon}$, \cref{eq:derivOmegaH} requires:
\begin{equation}\label{eq:hrocondi}
h_k =0,~ (k=1,2, \cdots, 6),
\end{equation}
in which $h_k$ are defined as:
\begin{subequations}
\label{eq:hyroCond}
\begin{align}
  \label{eq:hyroCond_a}
  &h_1=\int_{H_x} {\omega}^h \xi_1^2 \mathrm{d}V_{x'}-1, \\
  \label{eq:hyroCond_b}
  &h_2=\int_{H_x} {\omega}^h \xi_2^2 \mathrm{d}V_{x'}-1, \\
    \label{eq:hyroCond_c}
 &h_3=\int_{H_x} {\omega}^h \xi_3^2 \mathrm{d}V_{x'}-1, \\
    \label{eq:hyroCond_d}
  &h_4=\int_{H_x} {\omega}^h \xi_1\xi_2 \mathrm{d}V_{x'}, \\
    \label{eq:hyroCond_e}
  &h_5=\int_{H_x} {\omega}^h \xi_2\xi_3 \mathrm{d}V_{x'}, \\
    \label{eq:hyroCond_f}
  &h_6=\int_{H_x} {\omega}^h \xi_3\xi_1 \mathrm{d}V_{x'}.
\end{align}
\end{subequations}
To find the ${\omega}^h$ to satisfy \cref{eq:hrocondi} for arbitrary shapes of horizon $H_x$, we defined a function $I$ as:
\begin{equation}
    I(\omega^h,\lambda_1,\cdots, \lambda_6)=\frac{1}{2}\left({\omega^h-\omega^h_s} \right) \bullet \left({\omega^h-\omega^h_s} \right)- \sum_{k=1}^6 \lambda^h_k h_k,
\end{equation}
in which $\lambda^h_k$ are the Lagrange multipliers. The first variation of $I$ is:
\begin{equation}
\delta I = \left(\frac{\partial I}{\partial \omega^h } \right)^f \bullet \delta \omega^h + \sum_{k=1}^6 \frac{\partial I}{\partial \lambda^h_k} \delta \lambda^h_k,
\end{equation}
in which the superscript $f$ represents the Fréchet derivative and $\delta$ is the variation.
The condition $ \frac{\partial I}{\partial \lambda^h_k}=0$ leads to \cref{eq:hrocondi}.
The condition $\left(\frac{\partial I}{\partial \omega^h } \right)^f=0$ gives:
\begin{equation}\label{eq:omegah}
    \omega^h(\boldsymbol{\xi})=\omega^h_s(|\boldsymbol{\xi}|)+(\boldsymbol{\lambda}^h)^T \hat{\boldsymbol{\xi}}^h,
\end{equation}
in which $(\boldsymbol{\lambda}^h)^T \hat{\boldsymbol{\xi}}^h$ is the correction item for non-spherical horizon. $\boldsymbol{\lambda}^h$ and $\hat{\boldsymbol{\xi}}^h$ are defined as follows:
\begin{align}
    \boldsymbol{\lambda}^h & =\left[
    \begin{matrix}
    \lambda^h_1 & \lambda^h_2 & \lambda^h_3 & \lambda^h_4 & \lambda^h_5 & \lambda^h_6
    \end{matrix}
   \right]^T,\\
     \hat{\boldsymbol{\xi}}^h  & =\left[
    \begin{matrix}
    \xi_1^2 & \xi_2^2 & \xi_3^2 & \xi_1 \xi_2 & \xi_2 \xi_3 & \xi_3 \xi_1
    \end{matrix}
   \right]^T.
\end{align}
Substituting \cref{eq:omegah} into \cref{eq:hrocondi}, we have:
\begin{equation}\label{eq:lambdah}
   \boldsymbol{\lambda}^h=\left( \textbf{A}^h \right)^{-1}\textbf{R}^h,
\end{equation}
in which:
\begin{align}
    \textbf{A}^h & =\int_{H_x} \hat{\boldsymbol{\xi}}^h \otimes \hat{\boldsymbol{\xi}}^h \mathrm{d}V_{x'},\\
    \textbf{R}^h & =\textbf{r}^h- \int_{H_x} \omega_s^h \hat{\boldsymbol{\xi}}^h \mathrm{d}V_{x'},\\
    \textbf{r}^h & =\left[
    \begin{matrix}
    1 & 1 & 1 & 0 & 0 & 0
    \end{matrix}
 \right]^T.
\end{align}
The Lagrange multiplier $ \boldsymbol{\lambda}^h$ can be solved by \cref{eq:lambdah}. Then, the hydrostatic influence function is evaluated by \cref{eq:omegah}, which satisfies $\theta=\varepsilon_{ii}$ for arbitrary horizon shapes. It is worth noting that $\textbf{R}^h=0$, when $H_x$ is a sphere, leads to $\boldsymbol{\lambda}^h=0$ and $\omega^h=\omega^h_s$. Thus, the definition of the non-local dilatation of the XOSBPD in \cref{eq:theta_xpd} recovers that of the OSBPD.

Besides, for 2D problems, the vectors $\boldsymbol{\lambda}^h$, $\hat{\boldsymbol{\xi}}^h$ and $\textbf{r}^h$ are reduced to as:
\begin{align}
    \boldsymbol{\lambda}^h & =\left[
    \begin{matrix}
    \lambda^h_1 & \lambda^h_2 & \lambda^h_3
    \end{matrix}
   \right]^T,\\
    \hat{\boldsymbol{\xi}}^h & =\left[
    \begin{matrix}
    \xi_1^2 & \xi_2^2 & \xi_1 \xi_2
    \end{matrix}
   \right]^T,\\
    \textbf{r}^h & =\left[
    \begin{matrix}
    1 & 1 & 0
    \end{matrix}
 \right]^T.
\end{align}

\subsection{Deviatoric influence function}

Consider a 3-D body ($n_d=3$), we rearrange $ (\underline{\omega}^d \underline{e}^d )\bullet \underline{e}^d=\boldsymbol{\varepsilon}^d:\boldsymbol{\varepsilon}^d$ as:
\begin{equation}\label{eq:derivOmegaD1}
\begin{split}
    0&=(\underline{\omega}^d \underline{e}^d )\bullet \underline{e}^d - \boldsymbol{\varepsilon}^d:\boldsymbol{\varepsilon}^d\\
    &=\int_{H_x} {\omega}^d \underline{e}^d \underline{e}^d \mathrm{d}V_{x'} - \boldsymbol{\varepsilon}^d:\boldsymbol{\varepsilon}^d\\
    &=\int_{H_x} {\omega}^d \frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}^d\boldsymbol{\xi}}{|\boldsymbol{\xi}|} \frac{\boldsymbol{\xi}\cdot \boldsymbol{\varepsilon}^d\boldsymbol{\xi}}{|\boldsymbol{\xi}|} \mathrm{d}V_{x'} - \boldsymbol{\varepsilon}^d:\boldsymbol{\varepsilon}^d\\
    &=\varepsilon^d_{ij}\varepsilon^d_{mn}\int_{H_x} \frac{ \omega^d \xi_i \xi_j \xi_m \xi_n}{|\boldsymbol{\xi}|^2 } \mathrm{d}V_{x'} - \varepsilon^d_{ij}\varepsilon^d{ij}.
\end{split}
\end{equation}
Denote a 4th-order tensor $\mathbb{A}$ as:
\begin{equation}
    \mathbb{A}_{ijmn}=\int_{H_x} \frac{ \omega^d \xi_i \xi_j \xi_m \xi_n}{|\boldsymbol{\xi}|^2 } \mathrm{d}V_{x'},
\end{equation}
in which $\mathbb{A}$ is a symmetrical tensor:
\begin{equation}
    \mathbb{A}_{ijmn}=\mathbb{A}_{jimn}=\mathbb{A}_{ijnm}=\mathbb{A}_{jinm}=\mathbb{A}_{mnij}=\mathbb{A}_{nmij}=\mathbb{A}_{mnji}=\mathbb{A}_{nmji}.
\end{equation}
Using this symmetrical property and  $\varepsilon^d_{ij}=\varepsilon^d_{ji}$, \cref{eq:derivOmegaD1} can be expressed as:
\begin{equation}\label{eq:derivOmegaD2}
    \begin{split}
        0=&\varepsilon^d_{ij}\varepsilon^d_{mn} \mathbb{A}_{ijmn}- \varepsilon^d_{ij}\varepsilon^d_{ij} \\
        =&\left(\mathbb{A}_{1111} -1 \right) \varepsilon^d_{11}\varepsilon^d_{11}+\left(\mathbb{A}_{2222}-1\right) \varepsilon^d_{22}\varepsilon^d_{22}+\left(\mathbb{A}_{3333}-1\right)\varepsilon^d_{33}\varepsilon^d_{33} + \\
        & 2\mathbb{A}_{1122}\varepsilon^d_{11}\varepsilon^d_{22} + 2\mathbb{A}_{2233}\varepsilon^d_{22}\varepsilon^d_{33} + 2\mathbb{A}_{3311}\varepsilon^d_{33}\varepsilon^d_{11} + \\
        & \left(4\mathbb{A}_{1212}-2\right)\varepsilon^d_{12}\varepsilon^d_{12} + \left(4\mathbb{A}_{2323}-2\right)\varepsilon^d_{23}\varepsilon^d_{23} + \left(4\mathbb{A}_{1313}-2\right)\varepsilon^d_{13}\varepsilon^d_{13} + \\
      & 4\mathbb{A}_{1112}\varepsilon^d_{11}\varepsilon^d_{12} + 4\mathbb{A}_{1113}\varepsilon^d_{11}\varepsilon^d_{13} + 4\mathbb{A}_{2221}\varepsilon^d_{22}\varepsilon^d_{21} + 4\mathbb{A}_{2223}\varepsilon^d_{22}\varepsilon^d_{23} + 4\mathbb{A}_{3331}\varepsilon^d_{33}\varepsilon^d_{31} +
        4\mathbb{A}_{3332}\varepsilon^d_{33}\varepsilon^d_{32} + \\
    & 4\mathbb{A}_{1123}\varepsilon^d_{11}\varepsilon^d_{23} + 4\mathbb{A}_{2213}\varepsilon^d_{22}\varepsilon^d_{13} +
    4\mathbb{A}_{3312}\varepsilon^d_{33}\varepsilon^d_{12} + \\
    & 8\mathbb{A}_{1213}\varepsilon^d_{12}\varepsilon^d_{13} +                               8\mathbb{A}_{1223}\varepsilon^d_{12}\varepsilon^d_{23} +  8\mathbb{A}_{2313}\varepsilon^d_{23}\varepsilon^d_{13}.
    \end{split}
\end{equation}
At the first glance, for an arbitrary stain $\boldsymbol{\varepsilon}$, \cref{eq:derivOmegaD2} requires $\mathbb{A}_{1122}=0$ (see the 4th item) and $\mathbb{A}_{1212}=1/2$ (see the 7th item). However, from the symmetrical property of $\mathbb{A}$, we must have $\mathbb{A}_{1122}=\mathbb{A}_{1212}$. To solve this paradox, we first assume that:
\begin{equation}\label{eq:assum}
    \frac{1}{3}\mathbb{A}_{1111}=\frac{1}{3}\mathbb{A}_{2222}=\frac{1}{3}\mathbb{A}_{3333}=\mathbb{A}_{1122}=\mathbb{A}_{2233}=\mathbb{A}_{1133}.
\end{equation}
Using this assumption and $\varepsilon^d_{ii}=0$, we rearrange \cref{eq:derivOmegaD2} as:
\begin{equation}\label{eq:derivOmegaD3}
    \begin{split}
        0=&\left(\frac{2}{3}\mathbb{A}_{1111} -1 \right) \varepsilon^d_{11}\varepsilon^d_{11}+\left(\frac{2}{3}\mathbb{A}_{2222}-1\right) \varepsilon^d_{22}\varepsilon^d_{22}+\left(\frac{2}{3}\mathbb{A}_{3333}-1\right)\varepsilon^d_{33}\varepsilon^d_{33} + \\
        &\left(\frac{1}{3}\mathbb{A}_{1111} \varepsilon^d_{11}\varepsilon^d_{11} +\mathbb{A}_{1122} \varepsilon^d_{11}\varepsilon^d_{22} +\mathbb{A}_{1133} \varepsilon^d_{11}\varepsilon^d_{33} \right) +\\
         &\left(\mathbb{A}_{1122} \varepsilon^d_{11}\varepsilon^d_{22} +\frac{1}{3}\mathbb{A}_{2222} \varepsilon^d_{22}\varepsilon^d_{22} +\mathbb{A}_{2233} \varepsilon^d_{22}\varepsilon^d_{33} \right) +\\
         &\left(\mathbb{A}_{1133} \varepsilon^d_{11}\varepsilon^d_{33} +\mathbb{A}_{2233} \varepsilon^d_{22}\varepsilon^d_{33} +\frac{1}{3} \mathbb{A}_{3333} \varepsilon^d_{33}\varepsilon^d_{33} \right) +\\
        & \left(4\mathbb{A}_{1212}-2\right)\varepsilon^d_{12}\varepsilon^d_{12} + \left(4\mathbb{A}_{2323}-2\right)\varepsilon^d_{23}\varepsilon^d_{23} + \left(4\mathbb{A}_{1313}-2\right)\varepsilon^d_{13}\varepsilon^d_{13} + \\
      & 4\mathbb{A}_{1112}\varepsilon^d_{11}\varepsilon^d_{12} + 4\mathbb{A}_{1113}\varepsilon^d_{11}\varepsilon^d_{13} + 4\mathbb{A}_{2221}\varepsilon^d_{22}\varepsilon^d_{21} + 4\mathbb{A}_{2223}\varepsilon^d_{22}\varepsilon^d_{23} + 4\mathbb{A}_{3331}\varepsilon^d_{33}\varepsilon^d_{31} +
        4\mathbb{A}_{3332}\varepsilon^d_{33}\varepsilon^d_{32} + \\
    & 4\mathbb{A}_{1123}\varepsilon^d_{11}\varepsilon^d_{23} + 4\mathbb{A}_{2213}\varepsilon^d_{22}\varepsilon^d_{13} +
    4\mathbb{A}_{3312}\varepsilon^d_{33}\varepsilon^d_{12} + \\
    & 8\mathbb{A}_{1213}\varepsilon^d_{12}\varepsilon^d_{13} +                               8\mathbb{A}_{1223}\varepsilon^d_{12}\varepsilon^d_{23} +  8\mathbb{A}_{2313}\varepsilon^d_{23}\varepsilon^d_{13}\\
    =&\left(\frac{2}{3}\mathbb{A}_{1111} -1 \right) \varepsilon^d_{11}\varepsilon^d_{11}+\left(\frac{2}{3}\mathbb{A}_{2222}-1\right) \varepsilon^d_{22}\varepsilon^d_{22}+\left(\frac{2}{3}\mathbb{A}_{3333}-1\right)\varepsilon^d_{33}\varepsilon^d_{33} + \\
        & \left(4\mathbb{A}_{1212}-2\right)\varepsilon^d_{12}\varepsilon^d_{12} + \left(4\mathbb{A}_{2323}-2\right)\varepsilon^d_{23}\varepsilon^d_{23} + \left(4\mathbb{A}_{1313}-2\right)\varepsilon^d_{13}\varepsilon^d_{13} + \\
      & 4\mathbb{A}_{1112}\varepsilon^d_{11}\varepsilon^d_{12} + 4\mathbb{A}_{1113}\varepsilon^d_{11}\varepsilon^d_{13} + 4\mathbb{A}_{2221}\varepsilon^d_{22}\varepsilon^d_{21} + 4\mathbb{A}_{2223}\varepsilon^d_{22}\varepsilon^d_{23} + 4\mathbb{A}_{3331}\varepsilon^d_{33}\varepsilon^d_{31} +
        4\mathbb{A}_{3332}\varepsilon^d_{33}\varepsilon^d_{32} + \\
    & 4\mathbb{A}_{1123}\varepsilon^d_{11}\varepsilon^d_{23} + 4\mathbb{A}_{2213}\varepsilon^d_{22}\varepsilon^d_{13} +
    4\mathbb{A}_{3312}\varepsilon^d_{33}\varepsilon^d_{12} + \\
    & 8\mathbb{A}_{1213}\varepsilon^d_{12}\varepsilon^d_{13} +                               8\mathbb{A}_{1223}\varepsilon^d_{12}\varepsilon^d_{23} +  8\mathbb{A}_{2313}\varepsilon^d_{23}\varepsilon^d_{13}
    \end{split}
\end{equation}
For an arbitrary strain $\boldsymbol{\varepsilon}$, \cref{eq:derivOmegaD3} requires:
\begin{equation}\label{eq:devicondi}
    d_k=0, (k=1,2, \cdots, 15),
\end{equation}
in which $d_k$ is defined as:
\begin{align}\label{eq:odCondi}
\begin{matrix}
 d_1=\mathbb{A}_{1111}-\frac{3}{2},& d_2=\mathbb{A}_{2222}-\frac{3}{2},&  d_2=\mathbb{A}_{3333}-\frac{3}{2},\\
 d_4=\mathbb{A}_{1122}-\frac{1}{2},& d_5=\mathbb{A}_{2233}-\frac{1}{2},&  d_6=\mathbb{A}_{3311}-\frac{1}{2}, \\
 d_7=\mathbb{A}_{1112},& d_8=\mathbb{A}_{1113},& d_9=\mathbb{A}_{2221}, \\
d_{10}=\mathbb{A}_{2223},& d_{11}=\mathbb{A}_{3331},& d_{12}=\mathbb{A}_{3332}, \\
d_{13}=\mathbb{A}_{1123},& d_{14}=\mathbb{A}_{2213},&  d_{15}=\mathbb{A}_{3312}.
\end{matrix}
\end{align}
Note that \cref{eq:devicondi} satisfies and supports the assumption of \cref{eq:assum}. To find the deviatoric influence function $\omega^d$ that satisfies \cref{eq:devicondi} for arbitrary horizon shapes, we defined a function $J$ as:
\begin{equation}
    J(\omega^d,\lambda_1,\cdots, \lambda_{15})=\frac{1}{2}\left({\omega^d-\omega^d_s} \right) \bullet \left({\omega^h-\omega^h_s} \right)- \sum_{k=1}^{15} \lambda^d_k d_k,
\end{equation}
in which $\lambda^d_k$ are the Lagrange multipliers. The first variation of $J$ is:
\begin{equation}
\delta J = \left(\frac{\partial J}{\partial \omega^d } \right)^f \bullet \delta \omega^d + \sum_{k=1}^{15} \frac{\partial J}{\partial \lambda^d_k} \delta \lambda^d_k.
\end{equation}
The condition of $ \frac{\partial J}{\partial \lambda^d_k}=0$ leads to \cref{eq:devicondi}.
The condition of $\left(\frac{\partial J}{\partial \omega^d } \right)^f=0$ gives:
\begin{equation}\label{eq:omegad}
    \omega^d(\boldsymbol{\xi})=\omega^d_s(|\boldsymbol{\xi}|)+(\boldsymbol{\lambda}^d)^T \hat{\boldsymbol{\xi}}^d,
\end{equation}
in which $(\boldsymbol{\lambda}^d)^T\hat{\boldsymbol{\xi}}^d$ is the correction item for non-spherical horizons. $\boldsymbol{\lambda}^d$ and $\hat{\boldsymbol{\xi}}^d$ are defined as below:
\begin{equation}
\boldsymbol{\lambda}^d  =\left[
    \begin{matrix}
    \lambda^d_1 & \lambda^d_2 & \lambda^d_{3} & \cdots & \lambda^d_{15}
    \end{matrix}
   \right]^T,
\end{equation}
\begin{equation}
    \begin{split}
     \hat{\boldsymbol{\xi}}^d =\frac{1}{|\boldsymbol{\xi}|^2}&\left[
    \begin{matrix}
    \xi_1^4 & \xi_2^4 & \xi_3^4 & \xi_1^2 \xi_2^2& \xi_2^2 \xi_3^2 & \xi_3^2 \xi_1^2 & \xi_1^3 \xi_2 & \xi_1^3 \xi_3
    \end{matrix} \right.\\
    &
    \left.\begin{matrix}
    \xi_2^3 \xi_1 & \xi_2^3 \xi_3 & \xi_3^3 \xi_1 & \xi_3^3 \xi_2 & \xi_1^2 \xi_2 \xi_3 & \xi_2^2 \xi_1 \xi_3 & \xi_3^2 \xi_1 \xi_2
    \end{matrix}
   \right]^T.
\end{split}
\end{equation}
Substituting \cref{eq:omegad} into \cref{eq:devicondi}, we have:
\begin{equation}\label{eq:lambdad}
   \boldsymbol{\lambda}^d=\left( \textbf{A}^d \right)^{-1}\textbf{R}^d,
\end{equation}
in which:
\begin{align}
    \textbf{A}^d & =\int_{H_x} \hat{\boldsymbol{\xi}}^d \otimes \hat{\boldsymbol{\xi}}^d \mathrm{d}V_{x'},\\
    \textbf{R}^d & =\textbf{r}^d - \int_{H_x} \omega_s^d \hat{\boldsymbol{\xi}}^d \mathrm{d}V_{x'},\\
\textbf{r}^d & =\left[
    \begin{matrix}
    \frac{3}{2} & \frac{3}{2} & \frac{3}{2} & \frac{1}{2} & \frac{1}{2} & \frac{1}{2} & 0 & 0 & \cdots & 0
    \end{matrix}
 \right]^T_{15 \times 1}.
\end{align}
The Lagrange multiplier $ \boldsymbol{\lambda}^d$ can be solved by \cref{eq:lambdad}. Then, the deviatoric influence function can be evaluated by \cref{eq:omegad}, which satisfies $ (\underline{\omega}^d \underline{e}^d )\bullet \underline{e}^d=\boldsymbol{\varepsilon}^d:\boldsymbol{\varepsilon}^d$ for arbitrary horizon shapes. It is worth noting that $\textbf{R}^d=0$, when $H_x$ is a sphere, we have  $\boldsymbol{\lambda}^d=0$ and $\omega^d=\omega^d_s$. Thus, the definition of the non-local SED of the XOSBPD in \cref{eq:sde_xpd} recovers that of the OSBPD.

Besides, for 2D problems, the vectors $\boldsymbol{\lambda}^d$, $\hat{\boldsymbol{\xi}}^d$ and $\textbf{r}^d$ are reduced to:
\begin{align}
     \boldsymbol{\lambda}^d & =\left[
    \begin{matrix}
    \lambda^d_1 & \lambda^d_2 & \lambda^d_{3} & \lambda^d_{4} & \lambda^d_{5}
    \end{matrix}
   \right]^T,\\
     \hat{\boldsymbol{\xi}}^d & = \frac{1}{|\boldsymbol{\xi}|^2}\left[
    \begin{matrix}
    \xi_1^4 & \xi_2^4 & \xi_1^2 \xi_2^2 & \xi_1^3 \xi_2 & \xi_2^3 \xi_1
    \end{matrix}
   \right]^T,\\
\textbf{r}^d & =\left[
    \begin{matrix}
    \frac{3}{2} & \frac{3}{2} & \frac{1}{2} & 0 & 0
    \end{matrix}
 \right]^T.
\end{align}

\subsection{Scalar force state}
With the derived $\omega^h$ and $\omega^d$, the non-local SED is equal to the local SED:
\begin{equation}
     W_{pd}=W_{CCM}=\frac{\kappa \theta^2}{2}+\mu(\underline{\omega}^d \underline{e}^d )\bullet \underline{e}^d,
\end{equation}
where $\mu$ is the shear modulus and $\kappa$ is the bulk modulus defined as:
\begin{equation}
\kappa=
    \begin{cases}
    \frac{E}{3(1-2\nu)}, &\text{3D}\\
    \frac{E}{2(1+\nu)(1-2\nu)}, &\text{2D plane strain}\\
    \frac{E}{2(1-\nu)}, &\text{2D plane stress}\\
    \end{cases},
\end{equation}
where $E$ is the Young's modules and $\nu$ is the Poisson's ratio.

The scalar force state $\underline{t}(\textbf{x}) \langle \boldsymbol{\xi} \rangle$ for either spherical or non-spherical horizons can be derived by using the Fréchet derivative of \cref{eq:frederi}. However, \cref{eq:frederi} requires the prerequisites (see Eqs. (90) and (95) in \cite{silling_peridynamic_2007}):
\begin{subequations}
\label{eq:prereq}
\begin{align}
    \label{eq:prereq1}
    &\left(\frac{\partial W_{pd}}{\partial \underline{e}^h} \right)^f \bullet \Delta \underline{e}^d =0,\\
    \label{eq:prereq2}
    &\left(\frac{\partial W_{pd}}{\partial \underline{e}^d} \right)^f \bullet \Delta \underline{e}^h=0.
\end{align}
\end{subequations}

Using \cref{eq:eh_xpd,eq:ed_xpd,eq:theta_xpd,eq:hrocondi,eq:devicondi} and $ \varepsilon^d_{ii} = 0$, the prerequisites of \cref{eq:prereq1,eq:prereq2} for either spherical or non-spherical horizon can be proved as:
\begin{equation}
\begin{split}
    \left(\frac{\partial W_{pd}}{\partial \underline{e}^h} \right)^f \bullet \Delta \underline{e}^d &=\kappa \theta \left(\frac{\partial \theta }{\partial \underline{e}^h}\right)^f \bullet \Delta \underline{e}^d\\
    &=\kappa \omega^h |\boldsymbol{\xi}| \theta \bullet \Delta \underline{e}^d \\
    &= \int_{H_x} \kappa \omega^h |\boldsymbol{\xi}| \theta \frac{\Delta\varepsilon^d_{ij}\xi_i\xi_j}{|\boldsymbol{\xi}|} \mathrm{d}V_{x'}\\
    &= \kappa   \theta \Delta \varepsilon^d_{ij} \int_{H_x} \omega^h \xi_i\xi_j \mathrm{d}V_{x'}\\
    &=\kappa   \theta \Delta \varepsilon^d_{ii}\\ &=0,
\end{split}
\end{equation}
and
\begin{equation}
\begin{split}
  \left(\frac{\partial W_{pd}}{\partial \underline{e}^d} \right)^f \bullet \Delta \underline{e}^h
  &=\mu  \left(\frac{\partial \left(\underline{\omega}^d\underline{e}^d\bullet \underline{e}^d\right)}{\partial \underline{e}^d}\right)^f  \bullet \Delta \underline{e}^h \\
  &=2\mu \omega^d \underline{e}^d \bullet \Delta \underline{e}^h\\
  &= \int_{H_x} 2\mu \omega^d \frac{\varepsilon^d_{ij}\xi_i\xi_j}{|\boldsymbol{\xi}|} \frac{\Delta\varepsilon^h_{mn}\xi_m\xi_n}{|\boldsymbol{\xi}|} \mathrm{d}V_{x'}\\
    &= 2\mu \varepsilon^d_{ij} \Delta\varepsilon^h_{mn}  \int_{H_x} \omega^d \frac{\xi_i\xi_j\xi_m\xi_n}{|\boldsymbol{\xi}|^2} \mathrm{d}V_{x'}\\
    &=2\mu \varepsilon^d_{ij} \Delta\varepsilon^h_{mn} \mathbb{A}_{ijmn}\\
    &= 2\mu \varepsilon^d_{ij} \frac{\Delta\theta}{n_d} \mathbb{A}_{ijmm}\\
    &= \frac{2\mu \Delta \theta}{n_d} \frac{n_d+2}{n_d} \varepsilon^d_{ii} \\& = 0,
\end{split}
\end{equation}
in which $\theta$ is defined by \cref{eq:theta_xpd} and $\underline{e}^d$ is given by \cref{eq:ed_xpd}.

Therefore, the scalar force state $\underline{t}(\textbf{x}) \langle \boldsymbol{\xi} \rangle$ for either spherical or non-spherical horizon can be written as follows:
\begin{equation}\label{eq:t_xpd}
\begin{split}
    \underline{t}(\textbf{x}) \langle \boldsymbol{\xi} \rangle
     & =\left(\frac{\partial W_{pd}}{\partial \underline{e}^h}\right)^f + \left(\frac{\partial W_{pd}}{\partial \underline{e}^d}\right)^f \\
     &=\kappa \theta \left(\frac{\partial \theta }{\partial \underline{e}^h}\right)^f + \mu  \left(\frac{\partial \left(\underline{\omega}^d\underline{e}^d\bullet \underline{e}^d\right)}{\partial \underline{e}^d}\right)^f\\
     &=\kappa \omega^h |\boldsymbol{\xi}| \theta + 2\mu \omega^d \underline{e}^d.
\end{split}
\end{equation}
For a spherical horizon, we have $\omega^h=\omega^h_s$, $\omega^d=\omega^d_s$, and \cref{eq:t_xpd} recovers \cref{eq:finscalfor_sph}. Thus, the presented XOSBPD recovers the OSBPD when horizon is a sphere.

It is worth noting that, for some horizon shapes, $\omega^h$ and $\omega^d$ may be negative. For example, if all bonds in the horizon have $\xi_1 \xi_2>0$, there must exist at least one bond with $\omega^h<0$ because of $h_4=0$. It was stated in \citep{mitchell_position-aware_2015} that negative $\omega^h$ and $\omega^d$ are acceptable and do not introduce instability because they may not imply imaginary wave speeds. However, the authors believe that a negative $\omega^h$ is unacceptable because it violates physics. For example, the scalar force state should be negative when a body is under pure hydrostatic pressure with $\theta<0$ and $\underline{e}^d=0$. But, in this case, the evaluation of the scalar force state by \cref{eq:t_xpd} is positive if $\omega^h<0$. Furthermore, the authors found that the negative $\omega^h$ will also bring instability in the numerical implementation. To avoid this issue, the authors suggest that if there exists one bond with $\omega^h<0$ in the family, then for all bonds within this family, $\omega^h_s$ and $\omega^d_s$ have to be used to calculate the scalar force state. Note that $\omega^h$ is used to evaluate the $\theta$ to ensure its accuracy. Although this treatment introduces accuracy loss, it is acceptable because only a very few horizons have negative $\omega^h$. For example, a rectangular plate has only four horizons (centered at the four corners) with negative $\omega^h$.

Besides, for only a handful of cases, the determinants of the symmetrical matrices $\textbf{A}^h$ and $\textbf{A}^d$ may vanish. Slightly changing the horizon sizes or using a different weighted function $\omega(|\boldsymbol{\xi}|)$ may solve this problem.

\section{Numerical implementation}\label{sec:NumImp}

The static and dynamic problems can be solved with non-uniform discretization with various horizon sizes without volume and surface corrections by the XOSBPD. First, the Lagrange multipliers $\boldsymbol{\lambda}^h$ and $\boldsymbol{\lambda}^d$ are evaluated by \cref{eq:lambdah} and \cref{eq:lambdad}, respectively. Second, the influence functions $\omega^h$ and $\omega^d$ are calculated using \cref{eq:omegah} and \cref{eq:omegad}, respectively. Then, the force state is computed from \cref{eq:t_xpd,eq:T}. Finally, we compute the displacement results by the equation of motion of \cref{eq:eqm}. In this section, we briefly introduce the numerical implementation of the XOSBPD.

\subsection{Discretization}

\begin{figure}
    \centering
    \includegraphics[width=4in]{figures/varHorizon.png}
    \caption{Numerical implementation of the XOSBPD with non-uniform discretization, non-spherical horizons and various horizon sizes. $\delta_{(i)}$ and $\delta_{j}$ are different the horizon sizes. The exact horizons are determined by \cref{eq:horizon}.}
    \label{fig:varHorizon}
\end{figure}

The XOSBPD enables the non-uniform discretization with various horizon sizes. As shown in \cref{fig:varHorizon}, the domain is discretized into non-uniform elements by mesh generators, i.e., ANSYS. The element center is the PD node $\textbf{x}_{(i)}$ associated with the element volume of $V_{(i)}$ or area of $A_{(i)}$. Note that two PD nodes $\textbf{x}_{(i)}$ and $\textbf{x}_{(j)}$ must be paired interaction for either horizon. Otherwise, it will bring ghost forces and violate the conservation of linear and angular momentum. This issue can be remedied by the following family determination:
\begin{equation}\label{eq:horizon}
    H_{\textbf{x}_{(i)}}=\left\{\textbf{x}_{(j)} \in H_{\textbf{x}_{(i)}}: \left|\textbf{x}_{(j)}-\textbf{x}_{(i)} \right| \leq \delta_{(i)} \cup \left|\textbf{x}_{(j)}-\textbf{x}_{(i)} \right| \leq \delta_{(j)} \right\},
\end{equation}
where $\delta_{(i)}=m\Delta_{(i)} $ is the horizon size, $m$ is a constant, and $\Delta_{(i)}$ is defined as:
\begin{equation}
    \Delta_{(i)}=
    \begin{cases}
    \sqrt{A_{(i)}}, \qquad \text{2-D}\\
    \sqrt[3]{V_{(i)}}, \qquad \text{3-D}
    \end{cases}.
\end{equation}
By the meshfree discretization, the one point integration \citep{silling_meshfree_2005} scheme is used for the integral of field variable $f(\textbf{x},\textbf{x}')$ as:
\begin{equation}
    \int_{H_x} f(\textbf{x},\textbf{x}') \mathrm{d}V_{x'}=\sum_{j} f(\textbf{x},\textbf{x}_{(j)})V_{(j)}.
\end{equation}
The dynamic fracture of the PD problems can be solved by the explicit scheme as follows:
\begin{subequations}
\begin{align}
 \dot{\textbf{u}}^{n+1} &= \Ddot{\textbf{u}}^n \Delta t +  \dot{\textbf{u}}^{n}, \\
    {\textbf{u}}^{n+1} &= \dot{\textbf{u}}^{n+1} \Delta t +  {\textbf{u}}^{n}.
\end{align}
\end{subequations}
where $\Delta t$ is the time step and $n$ represents the step number.
\subsection{Adaptive dynamic relaxation}
Although the PD is essentially in a dynamic form, it also can be used to solve static or quasi-static problems using the adaptive dynamic relaxation (ADR) method \citep{Underwood_Dynamic_relaxation_1983}. The ADR approach introduces fictitious damping and inertia terms as:
\begin{equation}\label{eq:eqmADR}
    \textbf{M}\Ddot{\textbf{U}} + c \textbf{M} \dot{\textbf{U}}=\textbf{F},
\end{equation}
where \textbf{M} is the fictitious diagonal mass matrix, $\textbf{U}$ is the global displacement vector, $\textbf{F}$ is the resultant force vector, and $c$ is fictitious damping coefficient. \cref{eq:eqmADR} can be solved by the central-difference explicit method:
\begin{subequations}
\label{eq:adr}
\begin{align}
\label{eq:adr_a}
 \dot{\textbf{U}}^{n+1/2} &=\frac{(2-c^n\Delta t)\dot{\textbf{U}}^{n-1/2} + 2\Delta t \textbf{M}^{-1}
 \textbf{F}^n}{2+c^n \Delta t}, \\
 \label{eq:adr_b}
    {\textbf{U}}^{n+1} &= \dot{\textbf{U}}^{n+1/2} \Delta t +  {\textbf{U}}^{n}.
\end{align}
\end{subequations}
where the time step incremental is generally specified as $\Delta t= 1$. The damping coefficient is given as:
\begin{equation}
c^{n}=2 \sqrt{\frac{\left(\mathbf{U}^{n}\right)^{T}{ }^{1} \mathbf{K}^{n} \mathbf{U}^{n}} {\left(\mathbf{U}^{n}\right)^{T} \mathbf{U}^{n}}},
\end{equation}
in which ${ }^{1} \mathbf{K}^{n}$ is a diagonal matrix defined as:
\begin{equation}
{ }^{1} K_{i i}^{n}=\frac{ -\left(F_{i}^{n}-F_{i}^{n-1}\right) }{M_{i i}\Delta t \dot{U}_{i}^{n-1 / 2}}, \qquad (\text{No summation for the dummy index }i),
\end{equation}
where the diagonal mass matrix is selected by:
\begin{equation}
M_{i i} \geq \frac{1}{4} \Delta t^{2} \sum_{j}\left|K_{i j}\right|, \qquad (\text{No summation for the dummy index }i),
\end{equation}
and $\textbf{K}$ is the global stiffness matrix. For the OSBPD and XOSBPD, the global stiffness matrix is not explicitly given. However, the evaluation of $\textbf{K}$ does not require good accuracy for ADR. In this study, we use the linearized bond-based PD (LBBPD) with the small displacement assumption to calculate the $\textbf{K}$ for determine the diagonal mass $M_{ii} $, because the stiffness matrix of LBBPD is very easy to construct.

\subsection{The global stiffness matrix of the LBBPD}

For the BBPD, the bond force is only determined by the bond but not by any other bonds:
\begin{equation}
\textbf{f}(\boldsymbol{\xi}, \boldsymbol{\eta})=CS \frac{\boldsymbol{\xi}+ \boldsymbol{\eta}}{|\boldsymbol{\xi}+ \boldsymbol{\eta}|},
\end{equation}
in which $C$ is the material constant defined as:
\begin{equation}
C=\begin{cases}
\frac{6E}{\pi \delta ^3 (1-\nu)}, \qquad &\text{2-D plane stress}\\
\frac{6E}{\pi \delta ^3 (1+\nu)(1-2\nu)}, \qquad &\text{2-D plane strain}\\
\frac{18\kappa}{\pi \delta^4}, \qquad &\text{3-D}
\end{cases},
\end{equation}
and $S$ is the stretch defined as:
\begin{equation}
S=\frac{|\boldsymbol{\xi}+ \boldsymbol{\eta}|-|\boldsymbol{\xi}|}{|\boldsymbol{\xi}|}.
\end{equation}
Based on the small deformation assumption, the bond force can be linearized as \cite{prakash_multi-threaded_2020,zhang_ansys_2022}:
\begin{equation}
\textbf{f}(\boldsymbol{\xi}, \boldsymbol{\eta})=\frac{C}{|\boldsymbol{\xi}|^3}\left(\boldsymbol{\xi} \otimes \boldsymbol{\xi} \right) \boldsymbol{\eta}.
\end{equation}
Thus, the internal force applied on the points $\textbf{x}'$ and $\textbf{x}$ due to the bond $\boldsymbol{\xi}$ is:
\begin{equation}
\left[
    \begin{matrix}
    F(\textbf{x})\\
     F(\textbf{x}')
    \end{matrix}
    \right]=\textbf{K}_{bond}
    \left[
    \begin{matrix}
   \textbf{u} (\textbf{x})\\
     \textbf{u} (\textbf{x}')
    \end{matrix}
    \right],
\end{equation}
in which $\textbf{K}_{bond}$ is expressed as (3-D):
\begin{equation}
    \textbf{K}_{bond}=\left[
    \begin{matrix}
    -k_{11} & -k_{12}& -k_{13} &k_{11} & k_{12}& k_{13}\\
    -k_{21} & -k_{22}& -k_{23} &k_{21} & k_{22}& k_{23}\\
    -k_{31} & -k_{32}& -k_{33} &k_{31} & k_{32}& k_{33}\\
    k_{11} & k_{12}& k_{13} &-k_{11} & -k_{12}& -k_{13} \\
    k_{21} & k_{22}& k_{23} &-k_{21} & -k_{22}& -k_{23} \\
    k_{31} & k_{32}& k_{33} &-k_{31} & -k_{32}& -k_{33}
    \end{matrix}
    \right],
\end{equation}
where
\begin{equation}
k_{ij}=\frac{C\xi_i \xi_j}{|\boldsymbol{\xi}|^3}\mathrm{d}V_x \mathrm{d}V_{x'}.
\end{equation}
The global stiffness matrix of the LBBPD is assembled from the bond stiffness matrix as:
\begin{equation}
    \textbf{K}=\sum_{\text{all bonds}} \textbf{K}_{bond}.
\end{equation}

\begin{figure}[htbp]
    \centering
         \subfloat[]{\label{fig:plateHole
    \includegraphics[width=3.2in]{figures/plateHole.png}}
    \subfloat[]{\label{fig:plateHole_mesh
    \includegraphics[width=2.3in]{figures/plateHole_mesh.png}}
    \caption{A 2-D square plate with a central hole. The length of the square is $L=1$ m, the radius of the circle is $r=0.1$ m, and the displacement loading is $u_0 = 5 \times 10^{-4}$ m: (a) Geometry and loading, (b) Non-uniform meshes (7,680  non-uniform quadrilateral elements).}
    \label{fig:plate_hole_gemo_mesh}
\end{figure}

\begin{figure}[htbp]
    \centering
         \subfloat[]{\label{fig:plateHole_UX
    \includegraphics[width=2.6in]{figures/plateHole_UX.png}}\ \
    \subfloat[]{\label{fig:plateHole_UY
    \includegraphics[width=2.6in]{figures/plateHole_UY.png}}
    \caption{The displacement results of the 2-D plate with a central hole solved by the XOSBPD with ADR (a) $u_x$, (b) $u_y$.}
    \label{fig:plateHoleU}
\end{figure}

\begin{figure}[ht!]
    \centering
         \subfloat[]{\label{fig:comPlateHoleUx
    \includegraphics[width=2.8in]{figures/comPlateHoleUx.png}}\ \
    \subfloat[]{\label{fig:comPlateHoleUy
    \includegraphics[width=2.8in]{figures/comPlateHoleUy.png}}
    \caption{The displacement results along the edge of the hole of the 2-D plate by FEM (solid blue line), XOSBPD (dotted red line), and OSBPD (green solid line): (a) $u_x$, (b) $u_y$.}
    \label{fig:CompPlateHoleU}
\end{figure}

\section{Numerical results}\label{sec:NumR}

In this section, three classic examples used in the literature are simulated and compared. The first example is a 2-D plate with a central hole under displacement loading. The second example is extruding the 2-D plate in the $z$-direction with a depth of 0.3 m. The first two examples are static problems solved by the ADR to validate the accuracy of the presented XOSBPD for both 2-D and 3-D problems. The third example is the Kalthoff-Winkler experiment, a classic study of the dynamical fracture. The Kalthoff-Winkler experiment's simulation presented in this section is used to demonstrate the XOSBPD's capability of complex dynamical fracture analysis. For all the three examples, the weighted function is specified as $\omega(|\boldsymbol{\xi}|)=1$, and the horizon is determined by $\delta_{(i)}=3.01\Delta_{(i)}$.

\subsection{A 2-D plate with a central hole}\label{sec:plateHole}

\begin{figure}[ht!]
    \centering
         \subfloat[]{\label{fig:blockHole_UX
    \includegraphics[width=2.in]{figures/blockHole_UX.png}}\
    \subfloat[]{\label{fig:blockHole_UY
    \includegraphics[width=2.in]{figures/blockHole_UY.png}}\
    \subfloat[]{\label{fig:blockHole_UZ
    \includegraphics[width=2.in]{figures/blockHole_UZ.png}}
    \caption{The displacement results of the 3-D block with a central hole solved by the XOSBPD with ADR (a) $u_x$, (b) $u_y$, and (c) $u_z$.}
    \label{fig:blockHoleU}
\end{figure}

\begin{figure}
    \centering
         \subfloat[]{\label{fig:comBlockHoleUx
    \includegraphics[width=2.in]{figures/comBlockHoleUx.png}}\
    \subfloat[]{\label{fig:comBlockHoleUy
    \includegraphics[width=2.in]{figures/comBlockHoleUy.png}}\
    \subfloat[]{\label{fig:comBlockHoleUz
    \includegraphics[width=2.in]{figures/comBlockHoleUz.png}}\
    \caption{The displacement results along the edge of the hole of the 3-D block at $z=0.15$m by FEM (solid blue line), XOSBPD (dotted red line), and OSBPD (green solid line): (a) $u_x$, (b) $u_y$, (c) $u_z$.}
    \label{fig:CompBlockHoleU}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=2.5in]{figures/KWplate.png}
    \caption{Geometry and loading condition of the Kalthoff-Winkler’s experiment, with red lines representing the pre-existing cracks. The 2-D plane is subject to an impact load with a speed of $v_0 = 16.5$ m/s between the two pre-existing cracks (the red lines).}
    \label{fig:KWplate}
\end{figure}

\begin{figure}[ht!]
    \centering
         \subfloat[]{\label{fig:pwPlate24us
    \includegraphics[width=2.in]{figures/pwPlate24us.png}}\
    \subfloat[]{\label{fig:pwPlate36us
    \includegraphics[width=2.in]{figures/pwPlate36us.png}}\
    \subfloat[]{\label{fig:pwPlate48us
    \includegraphics[width=2.in]{figures/pwPlate48us.png}}\\
      \subfloat[]{\label{fig:pwPlate58
    \includegraphics[width=2.in]{figures/pwPlate58.4us.png}}\
    \subfloat[]{\label{fig:pwPlate75us
    \includegraphics[width=2.in]{figures/pwPlate75us.png}}\
    \subfloat[]{\label{fig:pwPlate104us
    \includegraphics[width=2.in]{figures/pwPlate104us.png}}\\
    \caption{Cracks evolution and deformation shape (scaled 5 times) of the Kalthoff-Winkler's experiment sample by the XOSBPD: (a) $t=24$ $\mu$s, th sample  of $\beta$ = 67° (respect to the $x$ axis). (b) $t=36$ $\mu$s, the two new cracks are growing in a linear manner. (c) $t=48$ $\mu$s, new damage emerges close to the center of the right-hand surface. (d) $t=58.4$ $\mu$s, the new damage branches into two cracks. (e) $t=75$ $\mu$s,  the two new cracks keep growing. and (f) $t=104$ $\mu$s, the two pre-existing and the new two cracks meet each other and the whole plane breaks into five pieces.}
    \label{fig:pwPlateCrack}
\end{figure}

\begin{figure}[ht!]
    \centering
     \subfloat[]{\label{fig:wp1DEM
    \includegraphics[width=3.in]{figures/wp1DEM.png}}\  \
    \subfloat[]{\label{fig:wp1XFEM
    \includegraphics[width=3.in]{figures/wp1XFEM.png}}\
    \caption{Cracks propagation simulation of the Kalthoff-Winkler's experiment using (a) DEM \cite{kosteski_crack_2012} and (b) XFEM \cite{belytschko_dynamic_2003}.}
    \label{fig:pwPlateCrackLite}
\end{figure}
As shown in Fig. \ref{fig:plateHole}, a 2-D plane stress plate with a central hole is subject to a displacement loading of $ u_0=5\times 10^{-4}\ \mathrm{m}$. The dimensions of the plate are $L=1$ m and $r=0.1$ m. The material properties of the plate are Young's modulus $E=70\ \mathrm{GPa}$ and Poisson's ratio $\nu=0.33$. The whole model is discretized into 7,680 non-uniform quadrilateral elements, and thus with a total of 7,680 PD nodes, as illustrated in \cref{fig:plateHole_mesh}.

\cref{fig:plateHoleU} illustrates the displacement results of this example solved by the XOSBPD with the ADR technique. We also simulate this example using the original OSBPD without volume and surface corrections and compare the results with FEM. \cref{fig:CompPlateHoleU} shows the displacement comparison results by the three methods along the edge of the hole. The arc angle $\beta$ along the edge is defined in \cref{fig:plateHole}. As shown in \cref{fig:CompPlateHoleU}, for both the horizontal displacement $u_x$ and the vertical displacement $u_y$, the XOSBPD is very close to the FEM while the OSBPD has a larger error. The displacement results of this 2-D example confirm that the presented XOSBPD works for the non-spherical horizon and non-uniform mesh.

\subsection{A 3-D block with a central hole}\label{sec:blockHole}

After the accuracy demonstration of the 2-D case of the XOSBPD, We extrude the 2-D plate of the first example in the $z$-direction with a depth of 0.3 m and simulate the 3-D block with the same loading conditions and same material properties as the first example. The 3-D block is discretized into 173,944 hexahedron elements, with a total number of 173,944 PD nodes.

\cref{fig:plateHoleU} illustrates domain discretization and the displacement results of this example solved by the presented XOSBPD with ADR technique.

Similar to the first example, \cref{fig:CompBlockHoleU} shows the displacement results by the FEM, XOSBPD, and original OSBPD methods along the edge of the hole on the surface $z=0.15$ m. As we can see, the displacement results of XOSBPD are very close to that of FEM. But for OSBPD, the results have significant errors, especially for $u_y$ and $u_z$. The displacement results of this 3-D example further confirm that the presented XOSBPD works for the non-spherical horizon and non-uniform mesh of 3-D problems.

\subsection{The Kalthoff-Winkler’s experiment} \label{sec:KW}

This section compares the dynamical fracture simulation of the well-known Kalthoff-Winkler's experiment \citep{kalthoff_modes_nodate_2000}. As depicted in \cref{fig:KWplate}, a 2-D plane strain plate is imposed by impact load with a speed of $v_0=16.5$ m/s between the two pre-existing cracks. The other boundaries are free. The material properties are: the Young's modulus $E=190$ GPa, the mass density $\rho=8,000\ \mathrm{kg/m^3}$, the Poisson's ratio $\nu=0.25$, and the critical fracture energy release rate $G_0=222,170\ \mathrm{J/m^2}$. A uniform mesh discretizes the domain with a size of 1.25 mm. The time step is specified as $\Delta t=80$ ns. For this example, the bond stretch criterion is used. That is, if the bond stretch $S$ exceeds the critical value $S_0$, the bond will be irreversibly broken. The critical stretch for plane strain is defined as \cite{zaccariotto_luongo_sarego_galvanetto_2015}:
 \begin{equation}
     S_0=\sqrt{\frac{5\pi G_0}{12E\delta}}.
 \end{equation}

The simulation results are presented in \cref{fig:pwPlateCrack}. As we can see, the two pre-existing cracks start to grow from around $t=24\ \mathrm{\mu s}$ with the inclined angle of 67° (respect to the $x$ axis), which is very close to 70° from the experimental observation \citep{kalthoff_modes_nodate_2000}. At around time $t=48\ \mathrm{\mu s}$, new damage emerges close to the center of the right-hand surface. At time $t=58.4\ \mathrm{\mu s}$, the new crack branches into two cracks and finally coalesce with the two old cracks (see \cref{fig:pwPlate104us}). The cracks growth trigger time ($t=24\ \mathrm{\mu s}$) is also obtain by \citet{kosteski_crack_2012}  using discrete element method (DEM) and by \citet{belytschko_dynamic_2003} using extended FEM (XFEM), and they also observed the new crack and its braching, as shown in \cref{fig:pwPlateCrackLite}.

\section{Conclusion}\label{sec:conc}

This work extends the OSBPD model from the spherical horizon to non-spherical ones. The formulation is achieved by introducing the Lagrange multipliers to ensure the non-local dilatation and non-local SED equal the local dilatation and local SED, respectively. Since the XOSBPD works for arbitrary horizon shapes, volume and surface corrections are no longer needed. Moreover, the non-uniform discretization implementation with various horizon sizes is also made possible by the XOSBPD, which can save computation costs and are more conformal to complex geometries. The first two static examples solved by the ADR show that the XOSBPD has excellent accuracy for both 2-D and 3-D problems. The third example compares with the Kalthoff-Winkler experiment and confirms the XOSBPD's complex dynamical fracture analysis capability. The proposed method paves a road to the comprehensive study of static and dynamics of the failure mechanism of structure components and solid materials.

\bibliographystyle{elsarticle-num-names}


\section{Introduction}
\label{sec:intro}
In microstructural optimization, one aims to find the optimal topology, within a representative unit cell, that maximizes the desired property.  This has many applications in engineering; for example,  energy dissipation \cite{asadpoure2017topology}, fluid applications \cite{guest2006optimizing}, thermal applications \cite{zhou2008computational}, phononic applications \cite{sigmund2003systematic}, medical implants \cite{hsieh2021architected}, and so on. Further, with the advent of additive manufacturing, the fabrication of such microstructures is possible today.

In a typical microstructural optimization problem, one attempts to maximize a quantity of interest (such as bulk modulus), subject to a mass constraint (or equivalently, volume-fraction constraint). However, in many applications, the desired mass is not known a priori.  Therefore, instead of imposing an arbitrary mass constraint, we consider imposing constraints on other physical quantities. The main objective of this paper is to develop a framework for solving such generalized microstructural optimization problems.

Further, it has been observed \cite {sigmund1994materials} that the design of negative Poisson ratio (NPR) materials, using standard optimization techniques, can be particularly challenging. Either specialized methods need to be developed \cite{yin2001optimality} or heuristic parameters must be used \cite{xia2015design}. Here, we show that standard L-BFGS optimization can be used for the robust design of NPR materials.

The remainder of this paper is organized as follows. First, the critical concept of homogenization is briefly reviewed in Section \ref{subsec:homogenization}. Then, the current methods of  microstructural optimization are reviewed in Section \ref{subsec:litreview}, with an emphasis on the classic density-based formulation.  Then, in Section \ref{subsec:generalizedProblems}, a generalized microstructural optimization problem. To solve such a class of problems, we introduce a neural-network (NN) framework in Section \ref{sec:method_DesignRepresentationNN} where the density field is represented through the weights associated with the NN. This allows for automatic sensitivity computation which is essential for solving generalized problems.  In Section \ref{subsec:augLag} the methodology to solve such generalized problems is discussed followed by the proposed algorithm  in Section \ref{subsec:algorithm}. Numerical examples are presented in Section \ref{sec:experiments}. Finally, conclusions and open issues are discussed in Section \ref{sec:conclusion}.

\section{Background}
\label{sec:background}

\subsection{Homogenization}
\label{subsec:homogenization}
In microstructural design, a common hypothesis is that the microstructure is locally periodic, and there is scale separation. For example, Figure \ref{fig:UnitCell} illustrates locally periodic microstructures and two representative unit cells.
\begin{figure}[H]
	\begin{center}
    \includegraphics[scale=0.5]{Figures/Method/unitcell.png}%
		\caption{Locally periodic microstructures and unit cells.}
		\label{fig:UnitCell}
	\end{center}
\end{figure}

Given a unit cell (microstructure) in 2D, a forward problem is to find its homogenized elasticity $3 \times 3$ tensor  $\boldsymbol{C}^H$. The theory of homogenization is well developed, for example, see \cite{yang2020determination, hassani1998review}. A typical numerical strategy \cite{andreassen2014determine} for computing  $\boldsymbol{C}^H$ is to impose three independent periodic boundary conditions (in 2D), and solve the resulting  finite element problems; see Figure \ref{fig:Homogenization}.  The stresses and strains from the three problems are then used to compute  $\boldsymbol{C}^H$ as described in \cite{andreassen2014determine}.

\begin{figure}[H]
	\begin{center}
    \includegraphics[scale=0.75]{Figures/Method/homogenization.png}%
		\caption{Extracting the elasticity tensor.}
		\label{fig:Homogenization}
	\end{center}
\end{figure}
Next, various quantities of interest, namely,  bulk modulus ($K$), shear modulus ($G$) and Poisson ratio  ($\nu$) can be extracted from $\boldsymbol{C}^H$ as follows:

\begin{subequations}
	\begin{align}
	 K = (C^H_{1,1} + C^H_{2,2} +C^H_{1,2} +C^H_{2,1})/4
    \label{eq:BulkModulus}\\
     G = C^H_{3,3}
    \label{eq:ShearModulus}\\
   \nu = (C^H_{2,1} + C^H_{1,2})/(C^H_{1,1} +C^H_{2,2})
    \label{eq:PoissonRatio}
	\end{align}
\end{subequations}
\subsection{Microstructural Optimization Methods}
\label{subsec:litreview}

The inverse problem is to arrive at an optimal topology that  maximizes or minimizes one of these quantities of interest. There are several methods available today for optimizing microstructures; see, for example, \cite{osanov2016topology, gao2020topology, vogiatzis2017topology, SureshMicrostructural, kollmann2020deep, guo2020semi}.

In the  popular density-based methods,  one defines a pseudo-density  $\rho_e\in(0,1]$ over the underlying finite element mesh. Then the microstructural optimization problem, of say, maximizing the bulk modulus, subject to a mass constraint (we prefer here a mass constraint over a volume constraint since this generalizes more easily to multiple materials), may be posed as follows:

\begin{subequations}
	\begin{align}
	& \underset{\boldsymbol{\rho}}{\text{minimize}}
	& &-K(\boldsymbol{\rho})\label{eqnObj_microstr}\\
	& \text{subject to}
	& & \mathbf{K}(\boldsymbol{\rho})\boldsymbol{u}^i = \boldsymbol{f}^i, i = 1,2,3\label{eqn:GoverningEqn_microstr}\\
	& & & \sum_e \rho_e v_e \lambda_e = \hat{m}\label{eqn:masscons_microstr}\\
	& & & 0 < \rho_e \le 1 \label{eqn:densityConstraint}
	\end{align}
\end{subequations}
where  $\mathbf{K}$ is the  stiffness matrix, $\boldsymbol{u}^i$ and $\boldsymbol{f}^i$ are the displacement vector and the external force vector for the three problems in Figure ~\ref{fig:Homogenization}, $v_e$ is the volume of the finite element,  $\lambda_e$ is the physical density of the base material, and $\hat{m}$ is the mass constraint. In addition, the solid isotropic material with penalization (SIMP) penalization model is employed to  link the density variables to the base material
\begin{equation}
    E(\rho_e) = E_{min} + E \rho_e^p
    \label{eq:SIMPModel}
\end{equation}
The field can now be optimized using, for example, optimality criteria \cite{Bendsoe2003} or MMA \cite{svanberg1987MMA}, resulting in the desired microstructural topology; see \cite{xia2015design}, for example. Similar problems can be posed, for example, to minimize the Poisson ratio, subject to a mass constraint.

\section{Proposed Framework}
\label{sec:methodology}

\subsection{Generalized Microstructural Problems}
\label{subsec:generalizedProblems}

As stated earlier, in many applications, the  mass constraint $\hat{m}$ (see Equation \ref{eqn:masscons_microstr}) is not known a priori.  Instead, it may be more advantageous to impose constraints on physical quantities. In this paper, we consider such generalized microstructural optimization problems of the form:

\begin{subequations}
	\begin{align}
	& \underset{\boldsymbol{\rho}}{\text{minimize}}
	& &\phi( \boldsymbol{C}^H (\boldsymbol{\rho}))\label{eqnObj_microstr}\\
	& \text{subject to}
	& & \boldsymbol{K}(\boldsymbol{\rho})\boldsymbol{u}^i = \boldsymbol{f}^i, i = 1,2,3\label{eqn:GoverningEqn_microstr}\\
	& & & h_j(\boldsymbol{\rho}) = 0, j = 1,2,... \label{eqn:equalityConstraint_microstr}\\
	& & & 0 < \boldsymbol{\rho} \le 1 \label{eqn:densityConstraint_microstr}
	\end{align}
\end{subequations}
where the objective $\phi( \boldsymbol{C}^H)$ represents one of the following: the bulk modulus ($K$), shear modulus ($G$), Poisson ratio  ($\nu$), or mass ($m$), and $h_j$ are equality constraints involving these quantities.

For example, the problem of maximizing the bulk modulus subject to a Poisson ratio constraint may be posed as:
\begin{subequations}
	\begin{align}
	& \underset{\boldsymbol{\rho}}{\text{minimize}}
	& &-K(\boldsymbol{\rho})\label{eqnObj_microstr}\\
	& \text{subject to}
	& & \boldsymbol{K}(\boldsymbol{\rho})\boldsymbol{u}^i = \boldsymbol{f}^i, i = 1,2,3\label{eqn:GoverningEqn_microstr}\\
	& & & \nu(\boldsymbol{\rho})/\hat{\nu} -1 = 0 \label{eqn:volcons_microstr}\\
	& & & 0 < \boldsymbol{\rho} \le 1 \label{eqn:densityConstraint_microstr}
	\end{align}
\end{subequations}
where $\hat{\nu}$ is the desired Poisson ratio. Similarly, the problem of minimizing the mass  subject to constraints on the bulk modulus and Poisson ratio  may be posed as:
\begin{subequations}
	\begin{align}
	& \underset{\boldsymbol{\rho}}{\text{minimize}}
	& &m(\boldsymbol{\rho})\label{eqnObj_microstr}\\
	& \text{subject to}
	& & \boldsymbol{K}(\boldsymbol{\rho})\boldsymbol{u}^i = \boldsymbol{f}^i, i = 1,2,3\label{eqn:GoverningEqn_microstr}\\
	& & & K(\boldsymbol{\rho})/\hat{K} -1 = 0 \label{eqn:GoverningEqn_microstr}\\
	& & & \nu(\boldsymbol{\rho})/\hat{\nu} -1 = 0 \label{eqn:volcons_microstr}\\
	& & & 0 < \boldsymbol{\rho} \le 1 \label{eqn:densityConstraint_microstr}
	\end{align}
\end{subequations}
where $\hat{K}$ is the desired bulk modulus, and $\hat{\nu}$ is the desired Poisson ratio. In the remainder of this paper, we will consider solving such problems.

\subsection{Representing Density using Neural Networks}
\label{sec:method_DesignRepresentationNN}

One of the challenges in solving such generalized problems is computing the sensitivities of the objective and constraints. Manual derivation  can be cumbersome and error-prone, especially when the problem is recast within the context of augmented Lagrangian formulation (see Section \ref{subsec:augLag}). We, therefore, propose here a neural-network (NN) framework that supports automatic sensitivity computation for gradient-based optimization \cite{ChandrasekharAuTO2021}. The framework not only eliminates the burden of manual sensitivity calculations, it offers other computational advantages as discussed in the remainder of the paper.

In particular, we employ a simple fully-connected feed-forward neural network \cite{bishop2006pattern}. The input to the network are points ($x$, $y$) within the domain; the NN has a series of hidden layers associated with activation functions  \cite{Goodfellow2016_deepLearning}, \cite{lu2019dyingReLU}. The output is the density  $\rho$ at that point; see Figure \ref{fig:NNarchitecture}. Observe that the final layer is  a SoftMax activation function  \cite{bishop2006pattern} that ensures that the  density lies between 0 and 1. The density  will depend on the weights and bias, i.e., they serve as the design variables $\boldsymbol{w}$.

\begin{figure}[H]
	\begin{center}
    \includegraphics[scale=0.5]{Figures/Method/singleMaterialFramework.png}%
		\caption{Neural network architecture for representing the density field.}
		\label{fig:NNarchitecture}
	\end{center}
\end{figure}

Thus, we can now repose the generalized problem as:

\begin{subequations}
	\begin{align}
	& \underset{\boldsymbol{w}}{\text{minimize}}
	& &\phi( \boldsymbol{C}^H (\boldsymbol{w}))\label{eq:eqnObj_microstr}\\
	& \text{subject to}
	& & \boldsymbol{K}(\boldsymbol{w})\boldsymbol{u}^i = \boldsymbol{f}^i, i = 1,2,3\label{eqn:GoverningEqn_microstr}\\
	& & & h_j(\boldsymbol{w}) = 0, j = 1,2,...n \label{eq:volcons_microstr}
	\end{align}
	\label{eq:general_NN_eqn}
\end{subequations}
Observe that an explicit (bound) constraint on the density field is not needed since it is automatically satisfied by the Softmax function.

To generalize the above framework  to multiple materials, the output of the NN is increased to $(M+1)$ variables, where $M$ is the number of non-void materials; see   Figure \ref{fig:NNarchitectureMultiple}. No other change is needed in the framework, i.e., one can use exactly the same number of design variables $\boldsymbol{w}$. Further, due to the nature of Softmax function, the partition of unity condition:
\begin{equation}
    \sum\limits_{m=0}^M  \rho_m = 1
    \label{eq:PartitionOfUnity}
\end{equation}
is  automatically satisfied, i.e., the sum of all densities is guaranteed to be unity. Finally, the SIMP material model in Equation  \ref{eq:SIMPModel} is generalized to multiple materials  as follows:
\begin{equation}
    E(\rho) = E_{min} +\sum\limits_{m=0}^M  E_m\rho_m^p
    \label{eq:SIMPMultipleMaterialModel}
\end{equation}

\begin{figure}[H]
	\begin{center}
    \includegraphics[scale=0.55]{Figures/Method/multiMaterialFramework.png}%
		\caption{Neural network architecture for representing multiple materials.}
		\label{fig:NNarchitectureMultiple}
	\end{center}
\end{figure}

\subsection{Augmented Lagrangian}
\label{subsec:augLag}

To solve both the single and multi-material  problems, we will use  the augmented Lagrangian method \cite{nocedal2006numericalOptimization}, i.e., let
\begin{equation}
L(\boldsymbol{w}) = \phi + \sum\limits_{j=1,2 ...}^{n} {\alpha_j} h_j^2 + \sum\limits_{j=1,2 ...}^{n} {\mu_j} h_j
\label{eq:lossEqn}
\end{equation}
where the penalty parameters $\alpha_j$ and Lagrange multipliers $\mu_j$ are updated through iterations. First, starting with a small value for $\alpha_j$ and a zero value for $\mu_j$  the augmented Lagrangian is minimized using the popular L-BFGS method \cite{nocedal2006numericalOptimization}. Then,  the coefficients $\alpha$ and $\mu$ are updated as follows (see Section \ref{sec:experiments}):
\begin{subequations}
	\begin{align}
	 {\alpha_j}^{(k+1)} =  {\alpha_j}^{(k)}  + \Delta \alpha
    \label{eq:augLagAlphaUpdate}\\
     {\mu_j}^{(k+1)}  = {\mu_j}^{(k)}  + 2{\alpha_j}^{(k+1)}  h_j^{(k)}
    \label{eq:augLagMuUpdate}
	\end{align}
\end{subequations}
The optimization is repeated until convergence. For termination, we consider two quantities
\begin{equation}
\epsilon_{\phi} = \mid (\phi^{(k+1)} - \phi^{(k)})/\phi^{(k+1)}\mid
  \label{eq:phiError}
\end{equation}
and
\begin{equation}
\epsilon_{h} = \sum\limits_{j=1,2 ...}^{n}\mid h_j \mid
\label{eq:gError}
\end{equation}
The algorithm terminates when both quantities are less than a prescribed value (see Section \ref{sec:experiments} for details). The overall framework is illustrated in Figure \ref{fig:optimizationFramework}.

\begin{figure}[H]
	\begin{center}
    \includegraphics[scale=0.65]{Figures/Method/optimizationFramework.png}%
		\caption{The optimization framework.}
		\label{fig:optimizationFramework}
	\end{center}
\end{figure}

\subsection{Algorithm}
\label{subsec:algorithm}

The algorithm for the proposed framework is described in \cref{alg:MicroTOuNN}. First, the target design objective $\phi$  and the set of equality constraints $h_j$ are specified along with an initial topology. Note that defining an initial topology is equivalent to initializing the neural network. Next, the domain is discretized into finite elements and the element centers serve as input to the neural network, which returns the elemental densities $\bm{\rho}$ as discussed in \cref{sec:method_DesignRepresentationNN}.  Homogenization is performed which predicts the homogenized structural response from the elemental densities as per \cref{subsec:homogenization}. Based on the results from homogenization, the objective $\phi$ \cref{eq:eqnObj_microstr} and the constraints $h_j$ \cref{eq:volcons_microstr} are evaluated.  Next, the loss function $L$ is computed using \cref{eq:lossEqn}, and the gradient of the loss function with respect to the neural network weights ($\nabla L$) is computed through automatic differentiation. The NN is then trained using PyTorch's L-BFGS optimizer as described in \cref{subsec:augLag}, which requires the update of Lagrangian parameters $\alpha,\mu$. The training procedure continues until the convergence criteria  are met.

\begin{algorithm}[H]
	\caption{Microstructural TO}
	\label{alg:MicroTOuNN}
	\begin{algorithmic}[1]
		\Procedure{MicroTOuNN}{$\Omega_0$, $\phi$, \ldots} \Comment{\parbox[t]{.3\linewidth}{Inputs}}

		\State $\Omega^0 \rightarrow \Omega^0_h$ \Comment{Domain discretization} \label{alg:domainDiscretize}

		\State $\bm{x} = \{x_e,y_e\}_{e \in \Omega^0_h} $ \Comment{elem centers; NN input} \label{alg:elemCenterComp}

		\State  k = 0; $\alpha_i = \alpha_0$; $\mu_i = 0 \; \forall i$ \Comment{Initialization (\cref{subsec:augLag})}
		\label{alg:Initialization }

		\Repeat \Comment{Optimization (Training)}

		\State $NN(\bm{x} ; \bm{w}) \rightarrow \bm{\rho}$ \Comment{Fwd prop NN} \label{alg:fwdPropNN}

        \State $ \bm{\rho} \rightarrow \boldsymbol{C}^H$ \Comment{Homogenization (\cref{subsec:homogenization})} \label{alg:homogenization}

        \State $ \boldsymbol{C}^H \rightarrow \phi $\Comment{Compute design objective \cref{eq:eqnObj_microstr}} \label{alg:objective}

  		\State $\{\boldsymbol{C}^H, \bm{\rho} \} \rightarrow h_j $ \Comment{Compute  constraints \cref{eq:volcons_microstr}} \label{alg:massConstraint}

		\State $\{\phi, h, \alpha_i, \mu_i\} \rightarrow L$ \Comment{Loss from \Cref{eq:lossEqn}} \label{alg:lossCompute}

		\State $AD(L \leftarrow \bm{w}) \rightarrow \nabla L $ \Comment{Auto diff for sensitivity} \label{alg:autoDiff}

		\State $\bm{w}  +  \Delta \bm{w} (\nabla L ) \rightarrow \bm{w} $ \Comment{L-BFGS  step} \label{alg:LBFGSStep}

		\State $\text{k}++$

		\State Update $\alpha_i, \mu_i$ \Comment {Aug-Lag terms \cref{eq:augLagAlphaUpdate}, \cref{eq:augLagMuUpdate}} \label{alg:OptPenaltyUpdate}

		\Until{ $ \Delta \phi < \hat{\epsilon}_\phi$ and $ \Delta h < \hat{\epsilon}_h $  and $k < k_{max}$} \Comment{convg. criteria}

	    \State	\Return $\bm{w}$

		\EndProcedure
	\end{algorithmic}
\end{algorithm}

\section{Numerical examples}
\label{sec:experiments}

In this section, we conduct several numerical experiments to illustrate the proposed algorithm. The implementation is in Python, within the PyTorch environment.  All experiments were conducted on an Intel i9-11900K @ 3.50GHz GHz, equipped with 32 GB of RAM. The default parameters are listed in Table \ref{table:defaultParameters}. We use the popular SIMP continuation scheme \cite{Rojas2015Continuation} for the material model. For the default NN configuration of 5 layers and 30 neurons per layer, the number of design variables $\bm{w}$ (weight and bias) is 3872.

\begin{table} [h!]
	\begin{center}
		\begin{tabular}{  r | l   }
			Parameter & Description and default value \\ \hline
			$E$, $\nu$ & For base material: $E=1$, $\nu=0.3$; for void: $E=0.001$, $\nu=0.3$  \\
		    $\lambda$(mass density) & For base material: $\lambda = 1$;  for void: $\lambda = 0.001$  \\
			NN &  Neural network: 5 layers and 30 neurons per layer, with Swish functions  \\
			$\alpha$, $\mu$ & Lagrangian parameters: $\alpha_0 = 1,  \Delta \alpha = 5$ and $\mu_0 = 0$\\
			$p$ & SIMP continuation parameters: $p_0 = 2$, $\Delta p = 0.5$, $ p_{max} = 10$ \\
			$n_x$,$n_y$ & Mesh elements: $n_x = 60$, $n_y = 60$ \\
			$\hat{\epsilon}_\phi$, $\hat{\epsilon}_h$ & Convergence criteria for objective and constraints: $\hat{\epsilon}_\phi$ = $\hat{\epsilon}_h$ = 0.025  \\

		\end{tabular}
	\end{center}
	\caption{Default simulation parameters.}
	\label{table:defaultParameters}
\end{table}
The default initial topology is a square block as illustrated in Figure \ref{fig:initialTopology}.

\begin{figure}[H]
	\begin{center}
    \includegraphics[scale=0.5]{Figures/Experiments/InitialTopology.png}%
		\caption{Default initial topology.}
		\label{fig:initialTopology}
	\end{center}
\end{figure}

 Through the experiments, we investigate the following:
\begin{enumerate}
	\item $\textit{Validation}$: First, we consider classic mass constrained optimization of bulk modulus, shear modulus, and Poisson ratio, and compare some of the computed values against well-established theoretical results.

	\item $\textit{Convergence}$: The typical convergence of the algorithm is then illustrated.

	\item $\textit{Initial Topology}$: Next, we consider initial topologies different from the one in Figure \ref{fig:initialTopology} and optimize for maximum bulk modulus.

	\item $\textit{Impact of NN Configuration}$: We then vary the NN size and consider its impact on Poisson ratio minimization.

	\item $\textit{Impact of Mesh Size}$: We repeat the above experiment but vary the mesh size instead of the NN size.

	\item $\textit{Generalized Problems}$: We then consider several generalized microstructural problems.

    \item $\textit{Multi-material design}$: Finally, we present results for multiple materials.

\end{enumerate}

\subsection{Validation}
\label{sec:expts_validation}

First, we consider the classic problem of  bulk modulus maximization subject to a mass-fraction constraint of 0.3. The resulting topology, computed in 47 seconds, is illustrated in Figure \ref{fig:maxBulkModulus}. Similarly,  Figure \ref{fig:maxShearModulus} illustrates a  microstructure with maximal shear modulus for the same mass constraint, computed in 53 seconds. Finally,  Figure \ref{fig:minPoissonRatio} illustrates the microstructure when the Poisson ratio is minimized for the same mass constraint, computed in 55 seconds,

\begin{figure}[H]
\begin{subfigure}[b]{.31\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/Validation/KValidation.png}%
	\end{center}
	\caption{Bulk modulus.}
	\label{fig:maxBulkModulus}
		\end{subfigure}
	\begin{subfigure}[b]{.31\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/Validation/GValidation.png}%
		\caption {Shear modulus.}
		\label{fig:maxShearModulus}
	\end{center}
	\end{subfigure}
		\begin{subfigure}[b]{.31\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/Validation/nuValidation.png}%
		\caption{Poisson ratio.}
		\label{fig:minPoissonRatio}
	\end{center}
	\end{subfigure}
	\caption{Classic microstructural optimization problems.}
\end{figure}
Table \ref{table:maxK} compares the computed bulk modulus for various mass fractions, against the Hashin-Shtrikman upper bound.
\begin{table}[ht]
	\centering
    \begin{tabular}{||l|l|l||}
			\hline
			$\hat{m}$ & $K^*$ &  $K_{HS}$  \\ \hline
			0.1 & 0.020 & 0.0282   \\ \hline
			0.3 & 0.092 & 0.099  \\ \hline
			0.5 & 0.18 & 0.20  \\ \hline
			0.7 & 0.34 & 0.35  \\  \hline
    \end{tabular}
    \caption{Computed bulk modulus vs. HS upper bound.  }
    \label{table:maxK}
\end{table}

\begin{comment}
\begin{equation}
    \boldsymbol{C}^H = \begin{bmatrix}
0.018 & -0.011 & 0\\
-0.011 & 0.018 & 0\\
0 & 0 & 0.0011\\
\end{bmatrix}
\end{equation}
\end{comment}

\subsection{Convergence}

Recall that the error in the objective $\epsilon _ \phi$ and error in the constraints $\epsilon _ g$, defined in Equation \ref{eq:phiError} and Equation\ref{eq:gError} respectively, are computed at the end of every L-BFGS iteration. The algorithm terminates when  $	\epsilon _ \phi  < \hat{\epsilon}_\phi$ \emph {and} $ \epsilon _ g < \hat{\epsilon}_g$ specified in Table \ref{table:defaultParameters}. Figure \ref{fig:KError} illustrates the objective error $\epsilon _ \phi $ and the constraint error $\epsilon _ g$  for bulk modulus maximization, while Figure \ref{fig:nuError} illustrates the two errors for Poisson ratio minimization.

\begin{figure}[H]
  \begin{subfigure}[b]{.45\textwidth}
  	\begin{center}
    \includegraphics[scale=0.4]{Figures/Experiments/Convergence/errorHistory_K.png}%
		\caption{Bulk modulus.}
		\label{fig:KError}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.45\textwidth}
  	\begin{center}
    \includegraphics[scale=0.4]{Figures/Experiments/Convergence/errorHistory_nu.png}%
		\caption{Poisson ratio.}
		\label{fig:nuError}
	\end{center}
	\end{subfigure}
	\caption{Objective and constraint errors after every L-BFGS iteration.}
	\label{fig:ErrorPlots}
\end{figure}

In the experiments, we observed that 3 to 20 L-BFGS iterations are sufficient. Further, each L-BFGS iteration typically involves 2 to 10 inner iterations. The homogenized properties are recorded at the end of every inner iteration. Figure \ref{fig:KConvergence} illustrates these quantities for bulk modulus maximization, and Figure \ref{fig:nuConvergence} illustrates these for Poisson ratio minimization. The spikes observed in the two figures correspond to the start of a new L-BFGS iteration when the penalty parameters are updated.

\begin{figure}[H]
  \begin{subfigure}[b]{.45\textwidth}
  	\begin{center}
    \includegraphics[scale=0.425]{Figures/Experiments/Convergence/convergence_K.png}%
		\caption{Bulk modulus.}
		\label{fig:KConvergence}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.45\textwidth}
  	\begin{center}
    \includegraphics[scale=0.425]{Figures/Experiments/Convergence/convergence_nu.png}%
		\caption{Poisson ratio.}
		\label{fig:nuConvergence}
	\end{center}
	\end{subfigure}
	\caption{Various homogenized quantities after every FE operation.}
	\label{fig:ErrorPlots}
\end{figure}

\subsection{Impact of Initial Topology}
\label{sec:initialTopology}
Next, we  consider initial designs different from Figure \ref{fig:initialTopology} and study their impact on the optimal topology, with all other parameters are kept constant at default values. Specifically, we revisit the problem of maximizing the bulk modulus, subject to a mass-fraction constraint of 0.3.

For the initial topology of a circular hole in Figure \ref{fig:circularHole}, the final topology is illustrated in  Figure \ref{fig:circularHoleMaxK}.
\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/ImpactOfInitialTopology/circularHole.png}%
		\caption{Initial topology.}
		\label{fig:circularHole}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
     \includegraphics[scale=0.3]{Figures/Experiments/ImpactOfInitialTopology/circularHoleMaxK.png}%
		\caption {Final topology.}
		\label{fig:circularHoleMaxK}
	\end{center}
	\end{subfigure}
	\caption{Maximizing bulk modulus with mass constraint.}
	\end{center}
\end{figure}
Similarly, with the initial topology of four circular holes in Figure \ref{fig:circularHoles}, the final topology, for the same problem is illustrated in  Figure \ref{fig:circularHolesMaxK}, with $K^*=0.090$.

\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/ImpactOfInitialTopology/circularHoles.png}%
		\caption{Initial topology.}
		\label{fig:circularHoles}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
     \includegraphics[scale=0.3]{Figures/Experiments/ImpactOfInitialTopology/circularHolesMaxK.png}%
		\caption {Final topology.}
		\label{fig:circularHolesMaxK}
	\end{center}
	\end{subfigure}
	\caption{Maximizing bulk modulus with mass constraint.}
	\end{center}
\end{figure}

Although the designs are different, the final bulk moduli are comparable. Thus, one can potentially discover different designs by changing the initial design.

\subsection{Impact of NN Size}
\label{sec:NNSize}
We now vary the neural-network size and minimize the Poisson ratio, subject to a mass-fraction constraint of 0.3 (with default parameters).  Figure \ref{fig:ImpactOfNN} illustrates various designs obtained for different NN configurations. The final Poisson ratio $\nu^*$ for the three cases are -0.68, -0.62 and -0.84 respectively. The number of design variables  are 492, 1782 and 11,682 respectively. We observe that as the number of design variables is increased, the algorithm tends to generate more complex designs. However, we did not observe any pattern between the NN size and the final objective achieved.

\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/NNSize/nu5x10.png}%
		\caption{NN: $5 \times 10$.}
		\label{fig:circularHoles}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/NNSize/nu5x20.png}%
		\caption {NN: $5 \times 20$.}
		\label{fig:circularHolesMaxK}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/NNSize/nu8x40.png}%
		\caption {NN: $8 \times 40$.}
		\label{fig:circularHolesMaxK}
	\end{center}
	\end{subfigure}
	\caption{Minimizing Poisson ratio with mass constraint using different NNs.}
	\label{fig:ImpactOfNN}
	\end{center}
\end{figure}

\subsection{Impact of Mesh Size}
\label{sec:MeshSize}

Next, we vary the mesh size and minimize the Poisson ratio, subject to a mass constraint of 0.3 (with all other parameters at default).  Figure \ref{fig:ImpactOfMeshSize} illustrates  the designs obtained for different mesh sizes. The final Poisson ratio $\nu^*$ for the three cases are -0.56, -0.80 and -0.65 respectively. Once again, we did not observe any pattern between the mesh size and the final objective achieved.

\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MeshSize/nu40x40.png}%
		\caption{Mesh: $40 \times 40$.}
		\label{fig:circularHoles}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MeshSize/nu80x80.png}%
		\caption {Mesh: $80 \times 80$.}
		\label{fig:circularHolesMaxK}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MeshSize/nu100x100.png}%
		\caption {Mesh: $100 \times 100$.}
		\label{fig:circularHolesMaxK}
	\end{center}
	\end{subfigure}
	\caption{Minimizing Poisson ratio for different mesh sizes.}
	\label{fig:ImpactOfMeshSize}
	\end{center}
\end{figure}

\subsection{Generalized Problems}
\label{sec:bulkModulusConstraint}
To illustrate the benefits of the proposed framework, we first used the MATLAB code published in \cite{xia2015design} to minimize the Poisson ratio to a mass constraint of $\hat{m} =0.5$. The algorithm did not terminate; the final topology, after a maximum allowable 450 FE operations (15 seconds), is illustrated in Figure \ref{fig:classic}; it exhibits the following characteristics $\nu = -0.23$, $K =0.015$ and $G = 0.002$. Observe that the microstructure is weak in bulk and shear.

In the proposed framework, the bulk modulus was maximized subject to $\nu = -0.23$ and $\hat{m} = 0.5$, i.e., the constraints are consistent with the results obtained above. The resulting topology, after 420 FE operations (52 seconds), illustrated in Figure \ref{fig:maxKGeneralized} exhibits the following characteristics $m =0.5$, $\nu = -0.23$, $K =0.043$ and $G = 0.006$. Thus the proposed framework increases the bulk and shear moduli by a factor of $3X$. For approximately the same number of finite element operations, the proposed framework is slower due to the overhead of automatic differentiation.

Next, the shear modulus was maximized subject to $\nu = -0.23$ and $\hat{m} = 0.5$. The resulting topology after 560 FE operations (91 seconds),  illustrated in Figure \ref{fig:maxGGeneralized} exhibits the following characteristics $m =0.5$, $\nu = -0.23$, $K =0.023$ and $G = 0.029$, i.e., the bulk and shear moduli increased by a factor of $10X$.

\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.5]{Figures/Experiments/Comparison/classic.png}%
		\caption{Mass constrained.}
		\label{fig:classic}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.45]{Figures/Experiments/Comparison/maxKGeneralized.png}%
		\caption {Maximizing $K$.}
		\label{fig:maxKGeneralized}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.45]{Figures/Experiments/Comparison/maxGGeneralized.png}%
		\caption {Maximizing $G$.}
		\label{fig:maxGGeneralized}
	\end{center}
	\end{subfigure}
	\caption{Classic versus proposed framework.}
	\label{fig:ClassicVersusProposed}
	\end{center}
\end{figure}

One can also impose multiple physical constraints within the proposed framework. Consider the minimization of the mass, subject to bulk modulus and Poisson ratio constraints. Figure \ref{fig:MultipleConstraints} illustrates three  designs obtained for three different scenarios. The corresponding mass fractions are 0.33, 0.69, and 0.79.
\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MultipleConstraints/ConstraintSet1.png}
		\caption{$\hat{K} = 0.1$, $\hat{\nu} = 0.1$.}
		\label{fig:circularHoles}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MultipleConstraints/ConstraintSet2.png}
		\caption{$\hat{K} = 0.1$, $\hat{\nu} = -0.1$.}
		\label{fig:circularHolesMaxK}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MultipleConstraints/ConstraintSet3.png}
		\caption{$\hat{K} = 0.1$, $\hat{\nu} = -0.2$.}
		\label{fig:circularHolesMaxK}
	\end{center}
	\end{subfigure}
	\caption{Minimizing mass with bulk modulus and Poisson ratio constraints.}
	\label{fig:MultipleConstraints}
	\end{center}
\end{figure}
One can also target a specific elasticity matrix $\boldsymbol{\hat{C}}^H$. We consider  a specific example from \cite{yin2001optimality}:
\begin{equation}
\boldsymbol{\hat{C}}^H = \begin{bmatrix}
0.13 & -0.03 & 0\\
-0.03 &  0.13 & 0\\
0 & 0 & 0.015
\end{bmatrix}
\label{eq:targetC}
\end{equation}
Note that there may exist multiple solutions (with significantly different masses) with the same $\boldsymbol{C}^H$ \cite{vogiatzis2017topology}. Using our framework by imposing the target $\boldsymbol{C}^H$ and a target mass as constraints, one can obtain different designs; these are illustrated in Figure \ref{fig:targetCSolutions}a and \ref{fig:targetCSolutions}b, with mass target of 0.37 and 0.51, respectively. A solution, with a mass of 0.37 was reported  in \cite{yin2001optimality}; see Figure \ref{fig:targetCSolutions}c.

\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/TargetC/targetSolution1.png}
		\caption{m = 0.37}
		\label{fig:target1}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/TargetC/targetSolution2.png}
	\caption{m = 0.51}
		\label{fig:target1}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.55]{Figures/Experiments/TargetC/yinSolution.png}
		\caption{m = 0.37 \cite{yin2001optimality}}
		\label{fig:reference}
	\end{center}
	\end{subfigure}
	\caption{Designs obtained for a targeted $\boldsymbol{\hat{C}}^H$.}
	\label{fig:targetCSolutions}
	\end{center}
\end{figure}

\subsection{Sampling at Higher Resolution}
\label{sec:multiMaterial}

Once the optimization is completed, the global representation of the density field via the NN allows us to sample the field at a finer resolution and  extract a high-resolution topology at no additional cost. This is illustrated in Figure \ref{fig:Resampling}. Observe that  directly smoothing the coarse topology will retain the original topological features. However, sampling at a high resolution can recover  topological features captured by the NN, as can be observed in Figure \ref{fig:Resampling}.

\begin{figure}[H]
\begin{center}
	\begin{subfigure}[b]{.4\textwidth}
  	\begin{center}
    \includegraphics[scale=0.4]{Figures/Experiments/NuConstrained/MaxKnu-0.75.png}%
		\caption {Raw output.}
	\end{center}
	\end{subfigure}
		\begin{subfigure}[b]{.4\textwidth}
  	\begin{center}
    \includegraphics[scale=0.4]{Figures/Experiments/NuConstrained/MaxKnu-0.75PostProcessed.png}%
		\caption {Sampling at higher resolution.}
	\end{center}
	\end{subfigure}

	\caption{Sampling the neural network at a higher resolution.}
	\label{fig:Resampling}
	\end{center}
\end{figure}

\subsection{Multiple Materials}
\label{sec:multiMaterial}

Next consider maximization of $K$, subject to a mass constraint ($\hat{m} = 0.3$), with two, three and four materials. The material properties are summarized in Table \ref{table:materialProp}.

\begin{table}[ht]
	\centering
\begin{tabular}{||l|l|l|l|l||}
			\hline
			Material & Color code & $E$  & $\nu$ & $\lambda$ \\ \hline
			0 & White & 0.001 & 0.3 & 0  \\ \hline
			1 & Black & 1 & 0.3 & 1  \\ \hline
			2 & Red & 0.2 & 0.3  & 0.2 \\ \hline
			3 & Green & 0.3 & 0.3  & 0.25 \\  \hline
			3 & Blue & 0.4 & 0.3 & 0.3 \\  \hline
\end{tabular}
\caption{Material properties.  }
\label{table:materialProp}
\end{table}

The results are illustrated in Figure \ref{fig:MultipleMaterials} (compare against Figure \ref{fig:maxBulkModulus} for single material). The final bulk modulus $K^*$ for the three cases are 0.172, 0.225 and 0.240, respectively (compared to 0.091 for single material). Observe that  the bulk modulus increases as additional materials are introduced.

\begin{figure}[H]
\begin{center}
  \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MultiMaterial/2MaterialsMaxK.png}%
		\caption{Two materials.}
	\end{center}
	\end{subfigure}
    \begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MultiMaterial/3MaterialsMaxK.png}%
		\caption {Three materials.}
	\end{center}
	\end{subfigure}
	\begin{subfigure}[b]{.3\textwidth}
  	\begin{center}
    \includegraphics[scale=0.3]{Figures/Experiments/MultiMaterial/4MaterialsMaxK.png}%
		\caption {Four materials.}
	\end{center}
	\end{subfigure}
	\caption{Maximizing $K$ with mass constraint using multiple materials. }
	\label{fig:MultipleMaterials}
	\end{center}
\end{figure}

\section{Conclusions}
\label{sec:conclusion}
This paper presents a generalized neural-network-based framework for microstructural optimization where any of the microstructural quantities, namely, bulk, shear, Poisson ratio, volume or mass, can serve as the objective, while the remaining can be subject to constraints.  The error-prone task of sensitivity computation was avoided by exploiting NN's backward propagation. Further, for designing NPR materials, we avoided the use of specialized optimization techniques and heuristics; instead, standard L-BFGS optimization was used. The framework was demonstrated using several numerical experiments. Due to the overhead cost of automatic differentiation (AD), the framework was found to be slower than, say, the MATLAB implementation presented in \cite{xia2015design}. However, we believe the benefits of AD outweigh the computational costs.

There are several directions for future research: extension to multi-physics \cite{das2020multi}, 3D \cite{Andreassen2014}, geometric and material non-linearity \cite{wallin2020nonlinear}, inclusion of stress constraints \cite{collet2018topology}, multi-stable materials \cite{yang20201d}, and inclusion of manufacturing constraints \cite{du2018connecting}.

\section*{Acknowledgments}
\label{sec:acknowledgements}

The authors would like to thank the support of the  National Science Foundation through grant CMMI 1561899, and  the U. S. Office of Naval Research under PANTHER award number N00014-21-1-2916  through Dr. Timothy Bentley.

\section*{Compliance with ethical standards}
\label{sec:ethics}
The authors declare that they have no conflict of interest.

\section*{Replication of Results}
\label{sec:replic}
The Python code pertinent to this paper is available at \href{https://github.com/UW-ERSL/MicroTOuNN}{https://github.com/UW-ERSL/MicroTOuNN}.