# nlfem: A flexible 2d FEM Code for Nonlocal Convection–Diffusion and
Mechanics

Manuel Klar Universitaet Trier, D-54286 Trier, Germany; Email: klar@uni-
trier.de<EMAIL_ADDRESS><EMAIL_ADDRESS>Christian
Vollmann11footnotemark: 1 Volker Schulz11footnotemark: 1

Abstract. In this work we present the mathematical foundation of an assembly
code for finite element approximations of nonlocal models with compactly
supported, weakly singular kernels. We demonstrate the code on a nonlocal
diffusion model in various configurations and on a two-dimensional bond-based
peridynamics model. Further examples can already be found in [17]. The code
nlfem is published under the MIT License111For details see, e.g.,
https://opensource.org/licenses/MIT. and can be freely downloaded at
https://gitlab.uni-trier.de/pde-opt/nonlocal-models/nlfem.

Keywords. Nonlocal operators, finite element discretizations, Python.

##  1 Introduction

Over the last two decades nonlocal models attracted attention due to their
capability of circumventing limitations of classical, local models based on
differential equations. While the theory and application of differential
equations is well established, nonlocal models still bear fundamental
questions. The models emerge due to various applications like anomalous
diffusion [7, 8], peridynamics [29] or image processing [20, 25] and are
diverse in their mathematical nature [16, 5].

Their investigation is often accompanied with numerical experiments motivated
by these applications. In this work, we describe the discretization of
nonlocal operators of the general form

$-\mathcal{L}_{\delta}\mathbf{u}(\mathbf{x}):=~{}2\int_{\widetilde{\Omega}\cap
B_{\delta}(\mathbf{x})}\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})d\mathbf{y},$

where the kernel $\mathbf{C}_{\delta}$ accounts for a finite range of nonlocal
interactions determined by the horizon $\delta>0$. The interaction
neighborhood $B_{\delta}(\mathbf{x})$, i.e., the support of the kernel, is
often given by a suitable approximation of the Euclidean or infinity norm
ball.

The purpose of nlfem is to compute numerical solutions to related boundary
value problems which run at convenient speed for researchers. The
discretization of these problems is achieved by a finite element
approximation. The resulting variational framework comes at the price of a
second integration layer compared to the operator in its strong form, which
makes it more costly compared to classical differential operators. Further
challenges in the implementation arise due to finite interaction horizons and
singularities of the kernel.

The development of efficient codes is an active field of research and there
are various related implementations. The foundation for operators related to
the fractional Laplacian222Throughout this work we always refer to the
fractional Laplacian in the integral form. is given by boundary element
methods [28]. Based on these fundamentals, a Matlab implementation for a
finite element approximation of the two–dimensional fractional Laplacian with
infinite interaction is presented in [1]. Further advanced techniques to
efficiently implement the fractional Laplacian and resulting dense stiffness
matrices are developed in [2] and incorporated into the finite element code
PyNucleus333https://github.com/sandialabs/PyNucleus. The latter is a
recommendable alternative to our code for nonlocal diffusion problems with
interaction neighborhoods based on the Euclidean norm ball. Specifically, it
can handle large dense matrices which makes it well–suited for a large or even
infinite horizon $\delta$.

Alternatives to kernel truncations have been investigated in [3]. Apart from
this, various structure exploiting approaches, finite-difference schemes and
meshfree methods have been used to compute numerical solutions. In particular
meshfree methods currently dominate the practical application of nonlocal
models for mechanics [24, 6]. However, implementations of finite element
methods for peridynamics have already been integrated into engineering
software [32]. For a general overview we refer the reader to the comprehensive
review paper [11] on numerical methods for nonlocal problems.

Our code nlfem assembles nonlocal operators on triangular meshes based on
linear continuous (CG) or discontinuous Galerkin (DG) ansatz spaces in 2d and
also offers this functionality in 1d and 3d for some generic settings.

It allows the assembly of problems with given Neumann and Dirichlet boundary
data [14, 27, 19]. An important detail is that the null space of systems
corresponding to pure Neumann boundary data is exact up to machine precision,
which allows an efficient evaluation of a pseudoinverse.

Concerning the domain, nlfem covers a variety of different configurations. It
handles nonlocal interactions in nonconvex or even disconnected domains where
the intersection between the interaction neighborhood $B_{\delta}(\mathbf{x})$
and the domain can be disconnected. For example, this is of particular
interest in shape optimization with nonlocal operators [31], where the domain
is modified iteratively.

The kernel can be symmetric or nonsymmetric as well as scalar– or
matrix–valued. For the symmetric case our discretization of the weak form
guarantees the symmetry of the stiffness matrix up to machine precision. The
code can generically handle smooth kernels and it comes with quadrature rules
for fractional type kernels as they are found in [28, 1]. In addition to that,
the kernel can depend on finite element labels which opens the door to a
monolithic assembly of interface problems determined by spatially variable
kernels.

The nlfem code is most efficient for operators with interaction horizons,
which are comparable to the mesh size, i.e., $h\leq\delta\leq Ch$ for some
$C\geq 1$. This relation is often used in the nonlocal mechanics setting; for
example, in [4, 23], the choices $\delta=3h$ and $\delta=4h$ are advocated,
respectively. In this regime the systems can be considered to be rather
sparse. A careful consideration of the quadrature and interpolation errors can
allow smaller rations $\delta/h$ and increases the sparsity of the related
systems.

Our implementation is based on the extensive discussions on the errors
incurred by various approximations of the interaction neighborhood and
quadrature rules [17, 31], all of which are implemented here. We only discuss
a small selection of interaction neighborhoods in this paper, such as two
approximations of the Eucldiean norm ball, which provably do not deteriorate
the finite element interpolation error [17]. Furthermore, the interaction
neighborhood of a kernel is efficiently determined by a breadth–first
traversal of finite elements throughout the assembly process, which avoids
expensive preprocessing computations.

A fundamental advantage of finite element methods is that they can be
considered to be asymptotically compatible in the sense of [30]. Our code
reproduces this property for the cases $h\leq\delta$, $h\sim\delta\to 0$ if
the implemented interaction neighborhood does not induce geometric errors. For
example, this is the case for the implemented infinity norm ball.

For convenience the assembly is performed in multiple threads and the main
routine, which is written in C++, comes with a user–friendly Python interface.
When it comes to solving, we note that the stiffness matrix is returned in
Compressed Sparse Row (CSR) format. Therefore, the user can apply any sparse
solver accessible from Python and apply it to the stiffness matrix.

The remainder of this article is organized into two main Sections. First, in
Section 2 we give a precise formulation of the targeted problem class in 2.1
along with its finite element approximation in 2.2. In its following
subsections we then highlight discretization details that deserve special
attention in a nonlocal framework. Second, in Section 3 we present various
numerical examples, including diffusion and mechanics, and give a brief
scaling study.

##  2 Finite Element Approximation

###  2.1 Problem Formulation

Let $\widetilde{\Omega}\subset\mathbb{R}^{d}$ be a compact base domain, and
let $\Omega\subset\widetilde{\Omega}$ be open in $\widetilde{\Omega}$444More
precisely, there exists an open set $\mathcal{O}\subset\mathbb{R}^{d}$, such
that $\mathcal{O}\cap\widetilde{\Omega}=\Omega$.. We note that the case
$\Omega=\widetilde{\Omega}$ is allowable, and $\Omega$ is open in
$\mathbb{R}^{d}$ only if
$\Omega\subset\operatorname{int}(\widetilde{\Omega})$. We refer to $\Omega$ as
domain, although it is not necessarily open in $\mathbb{R}^{d}$. The
complement of $\Omega$ in $\widetilde{\Omega}$ is denoted by
$\Omega_{D}:=\widetilde{\Omega}\setminus\Omega$. It can play the role of a
nonlocal Dirichlet boundary in suitable settings. If
$\widetilde{\Omega}=\Omega$ the resulting nonlocal problem can be interpreted
as Neumann type problem. In Figure 2.1 we present an exemplary configuration.
Let $\delta>0$ be an interaction horizon. We make the following assumption
about the kernel function.

######  Assumption 1.

We assume that for the matrix–valued _kernel function_
$\Psi_{\delta}:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{{n}\times{n}}$
there exists an $s\in(0,1)$ such that

$\displaystyle\Psi_{\delta}(\mathbf{x},\mathbf{y})\,\mathds{1}_{B_{\delta}(\mathbf{x})}(\mathbf{y})\,|\mathbf{x}-\mathbf{y}|^{{d}+2s}$
(1)

is bounded.

This allows singularities at the origin and also includes smooth kernels such
as the constant kernel. We assign to the above kernel function the interaction
neighborhood $B_{\delta}(\mathbf{x})$, where

$\displaystyle
B_{\delta}(\mathbf{x}):=\\{\mathbf{y}\in\mathbb{R}^{d}~{}|~{}|\mathbf{x}-\mathbf{y}|_{\bullet}\leq\delta\\},$
(2)

for some norm $|\cdot|_{\bullet}$ in $\mathbb{R}^{d}$. We denote the truncated
kernel by

$\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y}):=\Psi_{\delta}(\mathbf{x},\mathbf{y})\mathds{1}_{B_{\delta}(\mathbf{x})}(\mathbf{y}).$
(3)

The linear nonlocal operator under consideration acting on a function
$\mathbf{u}\colon\mathbb{R}^{d}\to\mathbb{R}^{n}$ is then given by

$\displaystyle-\mathcal{L}_{\delta}\mathbf{u}(\mathbf{x}):=p.v.~{}2\int_{\widetilde{\Omega}}\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})d\mathbf{y},$
(4)

where “p.v.” denotes the Cauchy principal value555More precisely,
$p.v.~{}\int_{\widetilde{\Omega}}h(\mathbf{x},\mathbf{y})d\mathbf{y}:=\lim_{\epsilon\to
0^{+}}\int_{\widetilde{\Omega}\setminus
B_{\epsilon}(\mathbf{x})}h(\mathbf{x},\mathbf{y})d\mathbf{y}$.. By testing (4)
with $\mathbf{v}\colon\widetilde{\Omega}\to\mathbb{R}^{n}$ where
$\mathbf{v}=0$ on $\Omega_{D}$, we obtain the bilinear form

$\displaystyle A(\mathbf{u},\mathbf{v})$
$\displaystyle:=-\int_{\Omega}\mathbf{v}(\mathbf{x})^{\\!\top\\!\\!}\mathcal{L}_{\delta}\mathbf{u}(\mathbf{x})d\mathbf{x}$
(5)
$\displaystyle=2\int_{\Omega}\mathbf{v}(\mathbf{x})^{\\!\top\\!\\!}~{}p.v.\int_{\widetilde{\Omega}}\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})d\mathbf{y}d\mathbf{x}.$

With $\mathbf{v}=0$ on $\Omega_{D}$ and Fubini we observe that

$\displaystyle\int_{\Omega}\mathbf{v}(\mathbf{x})^{\\!\top\\!\\!}~{}p.v.\int_{\widetilde{\Omega}}\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})d\mathbf{y}d\mathbf{x}$
$\displaystyle=$
$\displaystyle\int_{\widetilde{\Omega}}\mathbf{v}(\mathbf{x})^{\\!\top\\!\\!}~{}p.v.\int_{\widetilde{\Omega}}\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})d\mathbf{y}d\mathbf{x}$
$\displaystyle=$
$\displaystyle-\int_{\widetilde{\Omega}}~{}p.v.\int_{\widetilde{\Omega}}\mathbf{v}(\mathbf{y})^{\\!\top\\!\\!}\left(\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})\right)d\mathbf{y}d\mathbf{x}.$

So that the bilinear form can be written as

$\displaystyle A($ $\displaystyle\mathbf{u},\mathbf{v})$ (6)
$\displaystyle=\int_{\widetilde{\Omega}}\int_{\widetilde{\Omega}}(\mathbf{v}(\mathbf{x})-\mathbf{v}(\mathbf{y}))^{\\!\top\\!\\!}\left(\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})\right)d\mathbf{y}d\mathbf{x}$
$\displaystyle=\int_{\widetilde{\Omega}}\int_{\widetilde{\Omega}}\mathds{1}_{B_{\delta}(\mathbf{x})}(\mathbf{y})(\mathbf{v}(\mathbf{x})-\mathbf{v}(\mathbf{y}))^{\\!\top\\!\\!}\left(\Psi_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\Psi_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})\right)d\mathbf{y}d\mathbf{x}.$

Note that we exploited the symmetry of the exact indicator function
$\mathds{1}_{B_{\delta}(\mathbf{x})}(\mathbf{y})$ to obtain this equality.

###  2.2 Finite Dimensional Approximation

Let $\mathcal{T}^{h}:=\\{\mathcal{E}_{k}\\}_{k=1}^{K}$ denote a subdivision of
$\widetilde{\Omega}=\Omega\cup\Omega_{D}$ into polyhedral finite elements with
nodes $\\{\mathbf{x}_{m}\\}_{m=1}^{M}$.

######  Assumption 2.

We assume that $\Omega$ and $\Omega_{D}$ can be exactly covered by the
subdivisions
$\mathcal{T}^{h}_{\Omega}=\\{\mathcal{E}_{k}\\}_{k=1}^{K_{\Omega}}$ and
$\mathcal{T}^{h}_{D}=\\{\mathcal{E}_{k}\\}_{k=K_{\Omega}+1}^{K}$ with
$\mathcal{T}^{h}=\mathcal{T}^{h}_{\Omega}\cup\mathcal{T}^{h}_{D}$,
respectively, where $\mathcal{T}^{h}_{D}$ is possibly empty. Since we assume
polyhedral elements, this implies that

$\displaystyle\overline{\Omega}=\bigcup_{k=1}^{K_{\Omega}}\overline{\mathcal{E}}_{k}\textrm{~{}~{}~{}and~{}~{}~{}}\overline{\Omega}_{D}=\bigcup_{k=K_{\Omega}+1}^{K}\overline{\mathcal{E}}_{k}$
(7)

are polyhedral domains.

In the case of $\Omega_{D}=\emptyset$ we have that $K_{\Omega}=K$ and
$\mathcal{T}^{h}=\mathcal{T}^{h}_{\Omega}$. For convenience of notation we
assume an ordering of the nodes such that
$\\{\mathbf{x}_{m}\\}_{m=1}^{M_{\Omega}}\subset\Omega$ and
$\\{\mathbf{x}_{m}\\}_{m=M_{\Omega}+1}^{M}\subset\overline{\Omega}_{D}$. This
assumption is not made in the implementation of nlfem.

(a)

(b)

Figure 2.1:  (a) The compact base domain $\widetilde{\Omega}$ contains the
domains $\Omega$ (blue) and $\Omega_{D}$ (red), where the latter is possibly
empty. The blue boundary belongs to $\Omega$ and the red boundary belongs to
$\Omega_{D}$. (b) The red elements belong to the nonlocal Dirichlet boundary
$\Omega_{D}$, while the blue elements belong to the domain $\Omega$. Note that
the vertices on the red lines do not belong to the interior of $\Omega$.

We implement piecewise–linear continuous and discontinuous ansatz functions
$\\{\phi_{j}\\}_{{j}=1}^{J}$, where ${J}=nM$ in case of continuous and
${J}={n}({d}+1)K$ in case of discontinuous ansatz functions. Again for
convenience we assume an ordering of the basis functions and define the
corresponding finite dimensional subspaces

$\displaystyle
V^{h}(\widetilde{\Omega},\mathbb{R}^{n}):=\operatorname{span}(\\{\phi_{j}\\}_{{j}=1}^{{J}}),\textrm{~{}~{}and~{}~{}}V^{h}_{c}(\widetilde{\Omega},\mathbb{R}^{n}):=\operatorname{span}(\\{\phi_{j}\\}_{{j}=1}^{{J}_{\Omega}}).$
(8)

In case of continuous ansatz functions the unknown coefficients correspond to
the nodes lying in the interior of $\Omega$ with respect to
$\widetilde{\Omega}$. See Figure 2.1 for an illustration. Now, the evaluation
of the bilinear form $A$ on $V^{h}(\widetilde{\Omega},\mathbb{R}^{n})\times
V_{c}^{h}(\widetilde{\Omega},\mathbb{R}^{n})$ can be written as sum over the
finite elements, i.e.,

$\displaystyle A(\phi_{j},\phi_{i})=$ (9)
$\displaystyle\sum_{k=1}^{K}\sum_{\ell=1}^{K}\left[\int_{\mathcal{E}_{k}}\int_{\mathcal{E}_{\ell}}(\phi_{i}(\mathbf{x})-\phi_{i}(\mathbf{y}))^{\\!\top\\!\\!}\left(\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\phi_{j}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\phi_{j}(\mathbf{y})\right)d\mathbf{y}d\mathbf{x}\right].$

Since the kernel $\mathbf{C}_{\delta}$ may exhibit a truncation on some pairs
$(\mathcal{E}_{k},\mathcal{E}_{\ell})$ we need an appropriate approximation of
its support. Apart from a modeling aspect, a numerical advantage of the kernel
truncation is that it produces sparse systems on sufficiently coarse grids.
However, this can only be exploited if the numerical evaluation of the
interaction neighborhood $B_{\delta}(\mathbf{x})$ does not deteriorate the
interpolation error of the finite element subspace. While the approximation of
the infinity or Manhattan norm balls does not introduce a geometric error, the
one for curved neighborhoods, like the Euclidean norm ball, does.

In the following, we describe two major examples of the implemented ball
approximations for the Euclidean ball

$B_{\delta}^{2}(\mathbf{x}):=\\{\mathbf{y}\in\mathbb{R}^{d}~{}|~{}|\mathbf{x}-\mathbf{y}|_{2}\leq\delta\\},$

which we call the _approxcaps_ and the _nocaps_ approximations. Both are based
on the given finite element mesh $\mathcal{T}^{h}$. By “cap” we mean the
circular segments that arise when a finite element triangle is only partially
covered by the Euclidean ball; see Figure 2.2. Among others, these ball
approximations are investigated in [17].

(a)

(b)

Figure 2.2:  If a finite element triangle is only partially covered by an
Euclidean ball the intersection contains circular caps. The _nocaps_ ball (a)
omits this cap, whereas the _approxcaps_ ball (b) retriangulates the whole
intersection.

######  Definition 2.1 (_nocaps_ Ball).

For $\mathbf{x}\in\widetilde{\Omega}$ the _nocaps_ ball approximation is
defined as the convex hull of the intersection of the boundary $\partial
B^{2}_{\delta}(\mathbf{x})$ of the Euclidean ball and the boundaries of the
elements, i.e.,

$B_{\delta}^{ncp}(\mathbf{x}):=\operatorname{conv}\left(\bigcup_{\mathcal{E}_{\ell}\in\mathcal{T}^{h}}\partial\mathcal{E}_{\ell}\cap\partial
B_{\delta}^{2}(\mathbf{x})\right),$

where the dependency on the mesh size $h$ is omitted in the notation.

For ${d}=2$, one can show that the area of the symmetric difference of
$B^{2}_{\delta}(\mathbf{x})$ and $B_{\delta}^{ncp}(\mathbf{x})$ is provably of
order $\mathcal{O}(h^{2})$ on a given mesh of size $h$. This is error
commensurate with respect to the interpolation error of a linear finite
element ansatz space; details for the latter two statements can be found in
[17]. The convex hull omits circular caps which appear in the intersection of
some elements $\mathcal{E}_{\ell}\cap B_{\delta}^{2}(\mathbf{x})$ and have a
significant size on coarse grids. Therefore, by adding additional points on
the center of possible caps the geometric error can be reduced even further
while still being in $\mathcal{O}(h^{2})$ for $d=2$; see again Figure 2.2. The
maximum number of caps for a single intersection $\mathcal{E}_{\ell}\cap
B_{\delta}^{2}(\mathbf{x})$ is three.

While the results in [17] are derived for a fixed horizon $\delta$, related
investigations for the local limit ($\delta\to 0$) with polygonal ball
approximations can be found in [26].

######  Definition 2.2 (_approxcaps_ Ball).

Let $\mathbf{x}\in\widetilde{\Omega}$. We denote the points on the cap center
of each nonempty intersection $\mathcal{E}_{\ell}\cap\partial
B_{\delta}^{2}(\mathbf{x})$ by $\mathbf{y}_{\ell}$. The _approxcaps_ ball is
then defined by

$\displaystyle B_{\delta}^{acp}$ $\displaystyle(\mathbf{x}):=$
$\displaystyle\operatorname{conv}\left(B_{\delta}^{ncp}(\mathbf{x})\right.\cup\left.\,\\{\mathbf{y}_{\ell}\,|\,\mathbf{y}_{\ell}\textrm{
cap center of }\mathcal{E}_{\ell}\cap\partial
B_{\delta}^{2}(\mathbf{x})\textrm{ for some
}\mathcal{E}_{\ell}\in\mathcal{T}^{h}\\}\right).$

Exact quadrature rules for circular caps can be found in [18]. Here however,
the quadrature points have to be computed during run time, as the rules depend
on the geometry of the cap for higher quadrature orders. Therefore, we do not
consider these exact rules in nlfem.

In addition to these approximate balls, we next also introduce the infinity
normball.

######  Definition 2.3 (Infinity Norm Ball).

For $\mathbf{x}\in\widetilde{\Omega}$ the infinity normball is defined by

$B_{\delta}^{\infty}(\mathbf{x}):=\\{\mathbf{y}\in\mathbb{R}^{d}~{}|~{}|\mathbf{x}-\mathbf{y}|_{\infty}\leq\delta\\}.$

Proofs for the convergence of the nonlocal Dirichlet problem to the classical
Dirichlet problem with corresponding scaling constants for various kernel
functions can be found in [31]. Since the infinity normball is implemented
exactly, it allows numerical tests of the expected asymptotic compatibility of
the discretization scheme [30].

We finally note that the indicator function based on any implemented
truncation may lack in symmetry. More precisely, there might exist
$\mathbf{x}$ and $\mathbf{y}$, for which

$\displaystyle\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\neq\mathds{1}_{B_{\delta}^{\\#}(\mathbf{y})}(\mathbf{x}),$

where $B_{\delta}^{\\#}(\mathbf{x})$, $\\#\in\\{ncp,acp,\infty\\}$, represents
one of the implemented truncations. This artefact stems from the ball
approximation itself in case of the _nocaps_ and _approxcaps_ ball, but can
also be caused by the quadrature; see Remark 2.5 below.

We define the integrand

$\Phi_{i{j}}(\mathbf{x},\mathbf{y}):=(\phi_{i}(\mathbf{x})-\phi_{i}(\mathbf{y}))^{\\!\top\\!\\!}(\Psi_{\delta}(\mathbf{x},\mathbf{y})\phi_{j}(\mathbf{x})-\Psi_{\delta}(\mathbf{y},\mathbf{x})\phi_{j}(\mathbf{y}))$
(10)

and, based on the ball approximation and (6), the approximate bilinear form

$\displaystyle A_{h}^{\\#}(\phi_{j},\phi_{i})$
$\displaystyle:=\int_{\widetilde{\Omega}}\int_{\widetilde{\Omega}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\Phi_{i{j}}(\mathbf{x},\mathbf{y})d\mathbf{y}d\mathbf{x}$
(11)
$\displaystyle=\sum_{k=1}^{K}\sum_{\ell=1}^{K}\int_{\mathcal{E}_{k}}\int_{\mathcal{E}_{\ell}}\mathds{1}_{B_{\delta}^{\\#}}(\mathbf{x},\mathbf{y})\Phi_{i{j}}(\mathbf{x},\mathbf{y})d\mathbf{y}d\mathbf{x}.$

By Fubini’s integration theorem, for sufficiently smooth ansatz functions
$\phi_{j}\in V^{h}(\widetilde{\Omega},\mathbb{R}^{n})$ and $\phi_{i}\in
V_{c}^{h}(\widetilde{\Omega},\mathbb{R}^{n})$, the approximate bilinear form
defined in (11) can be written as

$\displaystyle A_{h}^{\\#}(\phi_{j},\phi_{i})$
$\displaystyle=\int_{\widetilde{\Omega}}\int_{\widetilde{\Omega}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\Phi_{i{j}}(\mathbf{x},\mathbf{y})d\mathbf{y}d\mathbf{x}$
$\displaystyle=\int_{\Omega}~{}p.v.\int_{\widetilde{\Omega}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\phi_{i}(\mathbf{x})^{\\!\top\\!\\!}\left(\Psi_{\delta}(\mathbf{x},\mathbf{y})\phi_{j}(\mathbf{x})-\Psi_{\delta}(\mathbf{y},\mathbf{x})\phi_{j}(\mathbf{y})\right)d\mathbf{y}~{}d\mathbf{x}$
$\displaystyle~{}~{}~{}-\int_{\widetilde{\Omega}}~{}p.v.\int_{\Omega}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\phi_{i}(\mathbf{y})^{\\!\top\\!\\!}\left(\Psi_{\delta}(\mathbf{x},\mathbf{y})\phi_{j}(\mathbf{x})-\Psi_{\delta}(\mathbf{y},\mathbf{x})\phi_{j}(\mathbf{y})\right)d\mathbf{y}~{}d\mathbf{x}$
$\displaystyle=2\int_{\Omega}\phi_{i}(\mathbf{x})^{\\!\top\\!\\!}~{}p.v.\int_{\widetilde{\Omega}}\mathds{1}^{S}_{B_{\delta}^{\\#}}(\mathbf{x},\mathbf{y})\left(\Psi_{\delta}(\mathbf{x},\mathbf{y})\phi_{j}(\mathbf{x})-\Psi_{\delta}(\mathbf{y},\mathbf{x})\phi_{j}(\mathbf{y})\right)d\mathbf{y}~{}d\mathbf{x},$

where

$\mathds{1}^{S}_{B_{\delta}^{\\#}}(\mathbf{x},\mathbf{y}):=\frac{1}{2}\left(\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})+\mathds{1}_{B_{\delta}^{\\#}(\mathbf{y})}(\mathbf{x})\right).$
(12)

In view of (5), this shows that the approximate bilinear form $A_{h}^{\\#}$
can be interpreted as the discretization of the operator
$-\mathcal{L}_{\delta}$ based on the symmetrified approximate indicator
function $\mathds{1}^{S}_{B_{\delta}^{\\#}}(\mathbf{x},\mathbf{y})$ instead of
$\mathds{1}_{B_{\delta(\mathbf{x})}}(\mathbf{y})$ in the strong form. However,
defining the approximate bilinear form as in (11) guarantees the symmetry of
the stiffness matrix $(A_{h}^{\\#}(\phi_{j},\phi_{i}))_{i,j}$; see also Remark
2.5 below.

In the following we discuss the evaluation of the local contributions

$\displaystyle
A^{\\#}_{{k}{\ell}}(\phi_{j},\phi_{i}):=\int_{\mathcal{E}_{k}}\int_{\mathcal{E}_{\ell}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\Phi_{i{j}}(\mathbf{x},\mathbf{y})d\mathbf{y}~{}d\mathbf{x}$
(13)

to the $(i,j)$-th entry of the stiffness matrix.

###  2.3 Population of the Stiffness Matrix

The assembly algorithm iterates over all pairs of elements and then adds the
local contributions (13) to the stiffness matrix. We call a pair of elements
$(\mathcal{E}_{k},\mathcal{E}_{\ell})$ _intersecting_ if

$\overline{\mathcal{E}_{k}}\cap\overline{\mathcal{E}_{\ell}}\neq\emptyset.$

If a pair is not intersecting, we call it _disjoint_. For the evaluation of
the contributions it is only important how the kernel behaves on a fixed pair
of elements. If for example the kernel exhibits a singularity at the origin
and the pair is disjoint, the singularity does not occur. Also, if a pair of
elements $(\mathcal{E}_{k},\mathcal{E}_{\ell})$ fulfills that

$\displaystyle\mathcal{E}_{\ell}\subset B_{\delta}(\mathbf{x})\textrm{ for all
}\mathbf{x}\in\mathcal{E}_{k},$

the truncation does not come into play. We therefore distinguish two cases. In
the first case the kernel has a singularity with $s\leq 0.5$ or the pair of
elements is disjoint. In the second case the kernel has a singularity with
$s>0.5$ and the elements are intersecting.

####  2.3.1 Disjoint Pairs or Kernels with $s\leq 0.5$

If $s\leq 0.5$ in (1) or if the pair $(\mathcal{E}_{k},\mathcal{E}_{\ell})$ is
disjoint, the local contributions (13) can be factored out and computed
separately, i.e.,

$\displaystyle A^{\\#}_{{k}{\ell}}(\phi_{j},\phi_{i})=$
$\displaystyle\int_{\mathcal{E}_{k}}\int_{\mathcal{E}_{\ell}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\phi_{i}(\mathbf{x})^{\\!\top\\!\\!}\Psi_{\delta}(\mathbf{x},\mathbf{y})\phi_{j}(\mathbf{x})d\mathbf{y}~{}d\mathbf{x}$
(14)
$\displaystyle-\int_{\mathcal{E}_{k}}\int_{\mathcal{E}_{\ell}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\phi_{i}(\mathbf{x})^{\\!\top\\!\\!}\Psi_{\delta}(\mathbf{y},\mathbf{x})\phi_{j}(\mathbf{y})d\mathbf{y}~{}d\mathbf{x}$
$\displaystyle+\int_{\mathcal{E}_{k}}\int_{\mathcal{E}_{\ell}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\phi_{i}(\mathbf{y})^{\\!\top\\!\\!}\Psi_{\delta}(\mathbf{y},\mathbf{x})\phi_{j}(\mathbf{y})d\mathbf{y}~{}d\mathbf{x}$
$\displaystyle-\int_{\mathcal{E}_{k}}\int_{\mathcal{E}_{\ell}}\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})\phi_{i}(\mathbf{y})^{\\!\top\\!\\!}\Psi_{\delta}(\mathbf{x},\mathbf{y})\phi_{j}(\mathbf{x})d\mathbf{y}~{}d\mathbf{x}.$
$\displaystyle=:$ $\displaystyle
A^{1}_{{k}{\ell}}(\phi_{j},\phi_{i})+A^{2}_{{k}{\ell}}(\phi_{j},\phi_{i})$
$\displaystyle+A^{3}_{{k}{\ell}}(\phi_{j},\phi_{i})+A^{4}_{{k}{\ell}}(\phi_{j},\phi_{i}).$

Therefore, nlfem allows continuous and discontinuous ansatz functions if
$s\leq 0.5$. Note, that the fractional Sobolev spaces
$H^{s}(\widetilde{\Omega})$ contain functions with jump discontinuities for
$0\leq s\leq 0.5$.

The expression $A^{1}_{{k}{\ell}}(\phi_{j},\phi_{i})$ is nonzero only if the
element $\mathcal{E}_{k}$ lies in the support of $\phi_{i}$ and $\phi_{j}$.
Similarly, the contribution of $A^{3}_{{k}{\ell}}(\phi_{j},\phi_{i})$ is
linked to the element $\mathcal{E}_{\ell}$. The term
$A^{2}_{{k}{\ell}}(\phi_{j},\phi_{i})$ is nonzero only if $\phi_{i}$ has its
support on $\mathcal{E}_{k}$ and $\phi_{j}$ on $\mathcal{E}_{\ell}$, where the
converse holds for $A^{4}_{{k}{\ell}}(\phi_{j},\phi_{i})$. That way we derive
the indices of the basis functions corresponding to the pair
$(\mathcal{E}_{k},\mathcal{E}_{\ell})$ in the stiffness matrix. We note that
the same indices of basis functions occur again for the pair
$(\mathcal{E}_{\ell},\mathcal{E}_{k})$, but the contributions to the stiffness
matrix are not identical as the truncation
$\mathds{1}_{B_{\delta}^{\\#}(\mathbf{x})}(\mathbf{y})$ is not symmetric.
Ultimately, the contribution of the two pairs
$(\mathcal{E}_{k},\mathcal{E}_{\ell})$ and
$(\mathcal{E}_{\ell},\mathcal{E}_{k})$ together lead to the symmetrified
truncation (12).

####  2.3.2 Kernels with $s>0.5$ on Intersecting Pairs

If a pair of elements is intersecting and $s>0.5$ the separation of integrands
in (14) is not admissible. Therefore, nlfem is restricted to continuous ansatz
functions if $s>0.5$. As the pair is intersecting, the elements are either
vertex touching, edge touching or identical. Therefore the number of ansatz
functions to be considered in these cases is 5, 4 or 3, respectively. If we
denote the vertices by $j_{1},\dots,j_{\ell}$ for $\ell=5,4,3$ we obtain 25,
16, or 9 pairs of ansatz functions $(\phi_{j_{\nu}},\phi_{j_{\nu^{\prime}}})$
for $\nu,\nu^{\prime}=1,\dots,\ell$ which yield nonzero contributions to the
local stiffness matrix
$(A^{\\#}_{{k}{\ell}}(\phi_{j_{\nu}},\phi_{j_{\nu^{\prime}}}))_{1\leq\nu,\nu^{\prime}\leq\ell}$.

###  2.4 Quadrature

The quadrature rules need to work for kernels with truncations and
singularites. However, the implemented quadrature rules in nlfem do not take
account of both at the same time on a fixed pair of elements. We therefore
require the following assumption.

######  Assumption 3.

The quadrature rules for singular kernels assume that for all intersecting
pairs of elements $(\mathcal{E}_{k},\mathcal{E}_{\ell})$ it holds that
$\mathcal{E}_{\ell}\subset B_{\delta}(\mathbf{x})$ for all
$\mathbf{x}\in\mathcal{E}_{k}$.

######  Remark 2.4.

Specific quadrature rules for singular kernels are required for sufficiently
strong singularities only. Kernels like the peridynamic kernel (26) can be
integrated by simply technically avoiding zero-divisions and we ignore those
cases in our discussion of singular kernels.

Given Assumption 3 and the fact that disjoint pairs do not require a treatment
with regularizing integral transforms we can evaluate all contributions as
given in Algorithm 1.

Algorithm 1 Evaluate $A^{\\#}_{{k}{\ell}}$ in (13)

1:if $\mathcal{E}_{k},\mathcal{E}_{\ell}$ intersect and the kernel is singular
then

2: Evaluate $A^{\\#}_{{k}{\ell}}$ via quadrature for singular kernels;

3: see Section 2.4.2 below

4:else

5: Evaluate $A^{\\#}_{{k}{\ell}}$ via quadrature for kernel truncations;

6: see Section 2.4.1 below

7:end if

In any case, the quadrature is performed by pulling back the domain of
integration to a reference domain

$\displaystyle\widehat{\mathcal{E}}\times\widehat{\mathcal{E}},~{}~{}\textrm{
where
}~{}~{}\widehat{\mathcal{E}}:=\\{\widehat{\mathbf{x}}\in\mathbb{R}^{d}~{}|~{}\widehat{\mathbf{x}}\geq
0,\sum_{\iota=1}^{{d}}\widehat{\mathbf{x}}_{\iota}\leq 1\\}.$

The affine linear mapping
$\chi_{k}:\widehat{\mathcal{E}}\rightarrow\mathcal{E}$ from the reference to a
physical element allows to define the pullback

$\displaystyle\widehat{\Phi}_{k\ell,i{j}}(\widehat{\mathbf{x}},\widehat{\mathbf{y}}):=\Phi_{i{j}}(\chi_{k}(\widehat{\mathbf{x}}),\chi_{\ell}(\widehat{\mathbf{y}})),$

the inner integral

$\displaystyle\widehat{\mathcal{K}}_{k\ell,i{j}}^{\\#}(\widehat{\mathbf{x}}):=\int_{\widehat{\mathcal{E}}}\mathds{1}_{B_{\delta}^{\\#}}(\chi_{k}(\widehat{\mathbf{x}}),\chi_{\ell}(\widehat{\mathbf{y}}))\widehat{\Phi}_{k\ell,i{j}}(\widehat{\mathbf{x}},\widehat{\mathbf{y}})d\widehat{\mathbf{y}}$
(15)

and the local contribution to the $({i},j)$-th entry of the stiffness matrix

$\displaystyle\widehat{A}_{k\ell,i{j}}^{\\#}:=\int_{\widehat{\mathcal{E}}}\widehat{\mathcal{K}}_{k\ell,i{j}}^{\\#}(\widehat{\mathbf{x}})d\widehat{\mathbf{x}}$
(16)

on the reference element.

####  2.4.1 Quadrature for Kernel Truncations

For some pairs $(\mathcal{E}_{k},\mathcal{E}_{\ell})$ we find that
$\mathcal{E}_{\ell}$ is only partially covered by the interaction neighborhood
$B_{\delta}^{\\#}(\mathbf{x})$ for some $\mathbf{x}\in\mathcal{E}_{k}$, so
that $\mathcal{E}_{\ell}\cap
B_{\delta}^{\\#}(\mathbf{x})\subsetneq\mathcal{E}_{\ell}$. In this case, the
ball approximations $B_{\delta}^{ncp}(\mathbf{x})$,
$B_{\delta}^{acp}(\mathbf{x})$ as well as the ball
$B_{\delta}^{\infty}(\mathbf{x})$ require a retriangulation of the integration
domain $\mathcal{E}_{\ell}\cap B_{\delta}^{\\#}(\mathbf{x})$. We denote the
set of elements which result from such a retriangulation by

$\displaystyle\mathcal{T}_{h,\ell}^{\\#}(\mathbf{x}):=\\{\tilde{\mathcal{E}}_{\ell}\\}_{\ell=1}^{L_{\ell}},$

so that $\bigcup\mathcal{T}_{h,\ell}^{\\#}(\mathbf{x})=\mathcal{E}_{\ell}\cap
B_{\delta}^{\\#}(\mathbf{x})$. Let
$\\{\widehat{\mathbf{y}}_{q},d_{{\widehat{\mathbf{y}}}_{q}}\\}_{q=1}^{Q}$
denote a quadrature rule on the reference element $\widehat{\mathcal{E}}$, and
let $\widehat{\mathbf{x}}\in\widehat{\mathcal{E}}$ be some reference point.
Then the fully-discrete inner integral from (15) reads as

$\widehat{\mathcal{K}}_{k\ell,i{j}}^{\\#}(\widehat{\mathbf{x}})\approx\sum_{\tilde{\ell}=1}^{L_{\ell}}\sum_{q=1}^{Q}\widehat{\Phi}_{k\tilde{\ell},i{j}}(\widehat{\mathbf{x}},\widehat{\mathbf{y}}_{q})d_{{\widehat{\mathbf{y}}}_{q}}.$
(17)

In finite element implementations the function values of ansatz functions are
usually pre-computed at the quadrature points and stored. However, if a
retriangulation is necessary the physical coordinates of the quadrature points
$\chi_{\tilde{\ell}}(\widehat{\mathbf{y}}_{q})$ for $q=1,\dots,Q$ on some
element
$\tilde{\mathcal{E}}_{\tilde{\ell}}\in\mathcal{T}_{h,\ell}^{\\#}(\mathbf{x})$
are known at runtime only and the ansatz functions are evaluated at the
corresponding points.

Now let
$\\{\widehat{\mathbf{x}}_{p},d_{{\widehat{\mathbf{x}}}_{p}}\\}_{p=1}^{P}$ be a
quadrature rule for the reference element of the outer integral. Then with
(17) the discretized version of the local contribution to the $(i,j)$–th entry
of the stiffness matrix (16) is obtained by

$\displaystyle\widehat{A}_{k\ell,i{j}}^{\\#}\approx\sum_{p=1}^{P}\left(\sum_{\tilde{\ell}=1}^{L_{\ell}}\sum_{q=1}^{Q}\widehat{\Phi}_{k\tilde{\ell},i{j}}(\widehat{\mathbf{x}}_{p},\widehat{\mathbf{y}}_{q})d_{{\widehat{\mathbf{y}}}_{q}}\right)d_{{\widehat{\mathbf{x}}}_{p}}.$

####  2.4.2 Quadrature for Singularities

Some kernel functions exhibit a singularity at the origin, which lies in the
integration domain of
$\widehat{\mathcal{K}}^{\\#}_{k\ell,i{j}}(\widehat{\mathbf{x}})$ whenever the
pair $(\mathcal{E}_{k},\mathcal{E}_{\ell})$ intersects. We therefore require
regularizing integral transforms [28]. Assumption 3 allows to ignore possible
truncations in (16) and to simply evaluate

$\displaystyle\widehat{A}_{k\ell,i{j}}^{\\#}=\int_{\widehat{\mathcal{E}}}\int_{\widehat{\mathcal{E}}}\widehat{\Phi}_{k\ell,i{j}}(\widehat{\mathbf{x}},\widehat{\mathbf{y}})d\widehat{\mathbf{y}}d\widehat{\mathbf{x}}.$
(18)

If Assumption 3 is violated, then the support of the kernel is overestimated
by (18). Intersecting pairs $(\mathcal{E}_{\ell},\mathcal{E}_{k})$ can be
vertex touching, edge touching or identical. For each of those cases we apply
integral transforms, which again pull back subsets of the integration domain
$\widehat{\mathcal{E}}\times\widehat{\mathcal{E}}$ to the unit cube
$(0,1)^{4}$. The transformations are well established and applied for example
in the field of boundary element methods. Details can be found in [28, 1].
Note again, that in the case of singular kernels with $s>0.5$ the implemented
routines cannot evaluate the expression in (18) for discontinuous basis
functions.

######  Remark 2.5 (Asymmetry due to quadrature).

We have mentioned in Section 2.2 that the approximate indicator function based
on the _approxcaps_ and _nocaps_ balls may generally lack in symmetry, and
therefore would lead to nonsymmetric stiffness matrices if one used
representation (5) of the bilinear form. Also truncations invoked by the
Infinity ball $B_{\delta}^{\infty}(\mathbf{x})$, which can be implemented
exactly, can be nonsymmetric on sufficiently irregular grids. These two
observations holds true independent of the symmetry of the kernel function
$\Psi_{\delta}(\mathbf{x},\mathbf{y})$. However, it is a desirable feature
that the symmetry of the kernel, i.e., the self-adjointness of the operator,
is transported through the discretization process. With other words, for a
symmetric kernel we expect a symmetric stiffness matrix. This behaviour is
more intricate and related to the approximation of the interaction domain of a
finite element by the union of interaction neighborhoods of quadrature points,
see Figure 2.3. Thus, we use by default representation (11) instead, which as
stated above corresponds in strong form to the operator
$-\mathcal{L}_{\delta}$ based on the approximate indicator function (12).
Thereby, we guarantee the stiffness matrix to be symmetric up to machine
precision for any symmetric kernel function. Also, an important consequence of
the symmetrification is that the null space of the stiffness matrix related to
pure Neumann type problems contains the constant vectors up to machine
precision. This allows the application of projected Krylov subspace methods or
rank-1 corrections to efficiently evaluate a pseudoinverse of the stiffness
matrix.

Figure 2.3: The Figure shows the approximation of the interaction domain of a
single triangular element by interaction neighborhoods of three quadrature
points for the infinity normball $B_{\delta}^{\infty}$.

We summarize the different cases for the quadrature of kernels in Table 1. The
table distinguishes between different degrees of singularity $s$ in case of
$d=2$ as well as intersecting and disjoint pairs of elements. When $s>0$ we
apply regularizing integral transforms represented by $\clubsuit$ in the
table. This quadrature rule is applied for intersecting pairs and requires
Assumption 3. Whenever $s\leq 0.5$ we allow for discontinuous and continuous
ansatz functions. When $s>0.5$ only continuous ansatz functions are provided.
For more details we refer to [2]. If $-1<s\leq 0$ the singularity might be so
weak, that it suffices to simply avoid zero-divisions on intersecting pairs,
which is represented by $\diamondsuit$ in the table. However, the quadrature
rule $\clubsuit$ can also be applied. If $s\leq-1$ there is no singularity at
all and symbol $\heartsuit$ represents a standard quadrature rule. The rules,
$\heartsuit$ and $\diamondsuit$, respect the kernel truncation, i.e. they do
not require Assumption 3.

| $s\leq-1$ | $-1<s\leq 0$ | $0<s\leq 0.5$ | $0.5<s<1$
---|---|---|---|---
Intersecting | $\heartsuit$ |  $\diamondsuit$ or $\clubsuit$ | $\clubsuit$ | $\clubsuit$
Disjoint | $\heartsuit$ | $\heartsuit$ | $\heartsuit$ | $\heartsuit$
| CG and DG | CG

$\heartsuit$ |  Standard quadrature rule on $\widehat{\mathcal{E}}\times\widehat{\mathcal{E}}$
---|---
$\diamondsuit$ |  Quadrature rule on $\widehat{\mathcal{E}}\times\widehat{\mathcal{E}}$ that avoids $\widehat{\mathbf{x}}_{p}=\widehat{\mathbf{y}}_{q}$
$\clubsuit$ |  Quadrature rule on $(0,1)^{4}$ after regularizing integral transforms

Table 1: The table gives an overview of the different quadrature rules which
can be applied for dimension $d=2$ in nlfem.

###  2.5 Traversal of the interaction neighborhood

The local contributions $\widehat{A}_{k\ell,i{j}}^{\\#}$ to the stiffness
matrix, defined in (13), can be nonzero even for pairs of remote finite
elements $(\mathcal{E}_{k},\mathcal{E}_{\ell})$, which we then refer to as
interacting elements in the following. Unstructured meshes and the various
supports of kernel functions call for a flexible routine to identify
interacting elements. The identification can be accomplished if the assembly
follows a breadth–first search. In order to describe the algorithm we define
the adjacency graph $\mathbb{T}_{adj}:=\left(\mathcal{T}^{h},E_{adj}\right)$
of the finite element mesh with vertices $\mathcal{T}^{h}$ and edges

$E_{adj}:=\\{(\mathcal{E}_{k},\mathcal{E}_{\ell})\in\mathcal{T}^{h}\times\mathcal{T}^{h}\,|\,\overline{\mathcal{E}_{k}}\cap\overline{\mathcal{E}_{\ell}}\neq\emptyset\\}.$

The graph $\mathbb{T}_{adj}$ can be understood as the dual graph of the finite
element mesh, and its vertices are therefore given by the elements. We
additionally define the interaction graph
$\mathbb{T}_{S}:=\left(\mathcal{T}^{h},E_{S}\right)$ with vertices
$\mathcal{T}^{h}$ and edges

$\displaystyle
E_{S}:=\left\\{\left.(\mathcal{E}_{k},\mathcal{E}_{\ell})\in\mathcal{T}^{h}\times\mathcal{T}^{h}\,\right|\,\widehat{A}_{k\ell,i{j}}^{\\#}\neq
0~{}\text{for some}~{}\phi_{j},\phi_{i}\in
V^{h}(\widetilde{\Omega},\mathbb{R}^{n})\right\\}.$ (19)

It is clear that $\mathbb{T}_{adj}$ is a spanning subgraph of $\mathbb{T}_{S}$
for any positive interaction horizon $\delta$. The adjacency graph
$\mathbb{T}_{adj}$ can be computed and stored efficiently, while the
interaction graph $\mathbb{T}_{S}$ exhibits storage requirements which are
comparable to the full stiffness matrix. We also note that the breadth–first
search described below allows to naturally identify intersecting pairs
$(\mathcal{E}_{k},\mathcal{E}_{\ell})$, as they only occur in the first layer
of the traversal. This important built–in feature is used to identify the
intersection cases for element pairs mentioned in Section 2.3.

Furthermore, for a fixed element $\mathcal{E}_{k}$ let us denote the set of
all interacting elements by
$\mathcal{T}^{h}_{k}:=\\{\left.\mathcal{E}_{\ell}\in\mathcal{T}^{h}\right|(\mathcal{E}_{k},\mathcal{E}_{\ell})\in
E_{S}\\}$. We then define the subgraph
$\mathbb{T}_{S_{k}}:=\left(\mathcal{T}^{h}_{k},E_{S_{k}}\right)$ of
$\mathbb{T}_{S}$ with vertices $\mathcal{T}^{h}_{k}$ and edges given by

$\displaystyle
E_{S_{k}}:=\left\\{\left.(\mathcal{E}_{k},\mathcal{E}_{\ell})\in\mathcal{T}^{h}\times\mathcal{T}^{h}\,\right|\,(\mathcal{E}_{k},\mathcal{E}_{\ell})\in
E_{S}\right\\}.$

######  Assumption 4.

We assume without loss of generality that each of the graphs
$\mathbb{T}_{S_{k}}$, $\mathbb{T}_{adj}$, and $\mathbb{T}_{S}$ is connected.

It is clear that Assumption 4 can be violated if $\widetilde{\Omega}$ is not
connected or even if it is connected as depicted in Figure 2.4(a). As a
remedy, the implementation allows to mark elements which are not part of the
finite element discretization. That way we can add artificial vertices and
elements to the mesh. One option then is to connect all vertices at the
boundary of the connected components of $\widetilde{\Omega}$ with one
artificial vertex, which introduces an adjacency between those elements.
Another, straightforward option is to embed $\widetilde{\Omega}$ into a
bounded and convex hold-all domain $\widehat{\Omega}\subset\mathbb{R}^{n}$,
which guarantees that the graphs $\mathbb{T}_{S_{k}}$, $\mathbb{T}_{adj}$, and
$\mathbb{T}_{S}$ are always connected, see Figure 2.4(b).

(a)

(b)

Figure 2.4:  (a) The interaction neighborhood of an element might be
disconnected even if the domain is connected. (b) A convex hold all domain
$\widehat{\Omega}$ allows to account for possible interactions within
disconnected interaction sets.

In that sense, Assumption 4 does not cause any loss of generality. If
Assumption 4 holds, we can recover the subgraph $\mathbb{T}_{S_{k}}$ with a
truncated breadth–first traversal of $\mathbb{T}_{adj}$ starting in the root
node $\mathcal{E}_{k}$ as given in Algorithm 2. To that end the root node
$\mathcal{E}_{k}$ is appended to a queue $Q:=[\mathcal{E}_{k}]$. While $Q$ is
not empty, an element $\mathcal{E}_{\widetilde{\ell}}$ is read and removed
from the queue. Its immediate neighbour, say
$\mathcal{E}_{\ell}\in\mathcal{N}(\mathcal{E}_{\widetilde{\ell}})$, are
obtained from the adjacency graph and the integrals
$\widehat{A}_{k\ell,i{j}}^{\\#}$ are successively evaluated for all
$\mathcal{E}_{\ell}\in\mathcal{N}(\mathcal{E}_{\widetilde{\ell}})$. The
element $\mathcal{E}_{\ell}$ is added to the queue whenever the integral does
not vanish, and finally marked as visited. This procedure automatically
truncates the search to interacting elements for any connected interaction
neighborhood.

Algorithm 2 Traversal of interaction neighborhood

1:for $\mathcal{E}_{k}$ in $\mathcal{T}^{h}$ in parallel threads do

2: Append $\mathcal{E}_{k}$ to the queue $Q$

3: while $Q$ is not empty do

4: Read and remove $\mathcal{E}_{\widetilde{\ell}}$ from $Q$

5: for $\mathcal{E}_{\ell}$ in local neighborhood
$\mathcal{N}(\mathcal{E}_{\widetilde{\ell}})$ do

6: if $\mathcal{E}_{\ell}$ is not visited then

7: Evaluate $A^{\\#}_{{k}{\ell}}$ (Algorithm 1)

8: if $A^{\\#}_{{k}{\ell}}\neq 0$ append $\mathcal{E}_{\ell}$ to $Q$

9: Mark $\mathcal{E}_{\ell}$ as visited

10: end if

11: end for

12: end while

13:end for

##  3 Numerical Examples

We solve a truncated fractional–type steady–state diffusion equation, a linear
bond–based peridynamics equation [29, 33] and a steady–state diffusion problem
based on the infinity normball. In the first two examples we demonstrate the
convergence rate of the approximate solutions to a manufactured solution as
the mesh size $h\to 0$. In the latter example we demonstrate the asymptotic
compatibility of the discretization scheme for a vanishing horizon $\delta$.

To this end let $\Omega:=(0,0.5)^{2}\subset\mathbb{R}^{2}$. We define the
interaction domain of $\Omega$ by
$\Omega_{D}:=[-\delta,~{}0.5+\delta]^{2}~{}\setminus~{}\Omega,$ so that
$\widetilde{\Omega}=[-\delta,0.5+\delta]^{2}$ is compact. We then want to
solve the nonlocal Dirichlet problem

$\displaystyle\begin{cases}-\mathcal{L}_{\delta}\mathbf{u}=\mathbf{f}&\textrm{on
}\Omega,\\\ \mathbf{u}=\mathbf{g}&\textrm{on }\Omega_{D}.\\\ \end{cases}$ (20)

We define the function spaces

$\displaystyle V(\widetilde{\Omega},\mathbb{R}^{n})=\\{\mathbf{u}\in
L^{2}(\widetilde{\Omega},\mathbb{R}^{n}):\|\mathbf{u}\|_{V}<\infty\\},$
$\displaystyle V_{c}(\widetilde{\Omega},\mathbb{R}^{n})=\\{\mathbf{u}\in
V(\widetilde{\Omega},\mathbb{R}^{n}):\mathbf{u}=0\textrm{ in }\Omega_{D}\\},$

where

$\displaystyle\|\mathbf{u}\|_{V}^{2}=A(\mathbf{u},\mathbf{u})+\|\mathbf{u}\|_{L^{2}(\widetilde{\Omega})}^{2}.$

For given data $\mathbf{f}\in L^{2}(\Omega,\mathbb{R}^{n})$ and
$\mathbf{g}:=\mathbf{v}_{|\Omega_{D}}$, where $\mathbf{v}\in
V(\widetilde{\Omega},\mathbb{R}^{n})$, we call $\mathbf{u}\in
V(\widetilde{\Omega},\mathbb{R}^{n})$ the weak solution to problem (20), if

$\displaystyle A(\mathbf{u},\mathbf{v})$
$\displaystyle=(\mathbf{f},\mathbf{v})~{}~{}\textit{for
all}~{}~{}\mathbf{v}\in V_{c}(\widetilde{\Omega},\mathbb{R}^{n}),$ (21)
$\displaystyle\textit{and}~{}~{}\mathbf{u}$ $\displaystyle=\mathbf{g}\textrm{
in }\Omega_{D}.$

The wellposedness of problem (21) for various choices of kernel functions can
be found, e.g., in [11, 13, 16, 15, 22]. By exploiting the known values of
$\mathbf{u}$ on the Dirichlet domain, we can rewrite (21) as

$\displaystyle
A_{\Omega\Omega}(\mathbf{u},\mathbf{v})=(\mathbf{f},\mathbf{v})-A_{\Omega\Omega_{D}}(\mathbf{g},\mathbf{v}),$
(22)

where

$\displaystyle A_{\Omega\Omega}(\mathbf{u},\mathbf{v}):=$
$\displaystyle\int_{\Omega}\int_{\Omega}(\mathbf{v}(\mathbf{x})-\mathbf{v}(\mathbf{y}))^{\\!\top\\!\\!}\left(\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})-\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{u}(\mathbf{y})\right)d\mathbf{y}d\mathbf{x}$
$\displaystyle+2\int_{\Omega}\int_{\Omega_{D}}\mathbf{v}(\mathbf{x})^{\\!\top\\!\\!}\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y})\mathbf{u}(\mathbf{x})d\mathbf{y}d\mathbf{x}$

and

$\displaystyle
A_{\Omega\Omega_{D}}(\mathbf{g},\mathbf{v}):=-2\int_{\Omega}\int_{\Omega_{D}}\mathbf{v}(\mathbf{x})^{\\!\top\\!\\!}\mathbf{C}_{\delta}(\mathbf{y},\mathbf{x})\mathbf{g}(\mathbf{y})d\mathbf{y}d\mathbf{x}.$

In the stiffness matrix the splitting (22) can be naturally obtained by
separating the columns corresponding to the degrees of freedom from the
columns corresponding to the nodes on the boundary $\Omega_{D}$.

###  3.1 Truncated Fractional–type Diffusion

We define the scalar–valued translationally invariant and symmetric kernel
function
$\gamma^{s}_{\delta}:\mathbb{R}^{2}\times\mathbb{R}^{2}\rightarrow\mathbb{R}$
by

$\displaystyle\gamma^{s}_{\delta}(\mathbf{x},\mathbf{y}):=c_{s,\delta}~{}\frac{1}{|\mathbf{x}-\mathbf{y}|^{2+2s}},$
(23)

where

$\displaystyle c_{s,\delta}=\frac{2-2s}{\pi\delta^{2-2s}}~{}~{}\textrm{ and
}~{}~{}s\in(0,1).$ (24)

Then the scalar–valued truncated fractional–type diffusion operator reads as

$\displaystyle-\mathcal{L}^{s}_{\delta}u(\mathbf{x}):=p.v.\int_{B_{\delta}(\mathbf{x})}\gamma^{s}_{\delta}(\mathbf{x},\mathbf{y})(u(\mathbf{x})-u(\mathbf{y}))d\mathbf{y}.$
(25)

The wellposedness of problem (21) for this choice of kernel is studied in
[14]. The constant $c_{\delta}$ depends on $\delta$ and $s$ and is chosen such
that the operator converges to the classical Laplacian as $\delta\to 0$; see,
e.g., [31, Lemma 7.4.1]. For another choice of the constant convergence to the
fractional Laplacian as $\delta\to\infty$ can also be obtained [12].

In the example above, we choose the manufactured solution
$u(\mathbf{x})=\mathbf{x}_{1}^{2}\mathbf{x}_{2}+\mathbf{x}_{2}^{2}$ and set
$f(\mathbf{x}):=-\Delta u(\mathbf{x})=-2(\mathbf{x}_{2}+1)$ in $\Omega$ and
$g(\mathbf{x}):=u(\mathbf{x})$ on $\Omega_{D}$. Since the correctly scaled
nonlocal operator equals the classical Laplacian operator on polynomials of
order up to three (see, e.g., [31]), we have that $u(\mathbf{x})$ is the
solution of problem (21). Furthermore, we choose $s:=0.5$, $\delta=0.2$ and
various mesh sizes $h$ as given in the tables below. In view of Table 1, for
pairs of disjoint elements we use as $\heartsuit$ a 7-point quadrature
rule666Specifically, the quadrature points are the barycenter, the vertices,
and the mid-side points of the reference triangle $\widehat{\mathcal{E}}$ and
the corresponding weights are $\frac{27}{60}\\!\cdot\\!\frac{1}{2}$,
$\frac{3}{60}\\!\cdot\\!\frac{1}{2}$, and
$\frac{8}{60}\\!\cdot\\!\frac{1}{2}$, respectively. on each, i.e., outer and
inner, element. Since six of the seven points are located on the boundary of
the triangle, this choice has proven to be advantageous in case of truncated
kernels as it mitigates the issue described in Figure 2.3; for more details
see [17]. For intersecting pairs we need a quadrature rule on $(0,1)^{4}$
after the integral transformations are performed. For $\clubsuit$ we choose a
tensor product of a 5-point Gauss quadrature rule.

The convergence rates on a continuous Galerkin ansatz space are shown in Table
2. For our choice of $s=0.5$ a discontinuous Galerkin ansatz space is also
conforming, and the results are presented in Table 3. In both settings, we
observe the expected second–order convergence rate as the mesh size $h\to 0$;
see, e.g., [14]. Note that the given examples violate Assumption 4 in the
first stage of the experiments as $2h>\delta$. We see that the first rates in
both tables are affected by this.

$h$ | dof | L2 Error | Rates
---|---|---|---
1.41e-01 | 1.60e+01 | 7.19e-04 | 0.00e+00
7.07e-02 | 8.10e+01 | 1.65e-04 | 2.13e+00
3.54e-02 | 3.61e+02 | 4.09e-05 | 2.01e+00
1.77e-02 | 1.52e+03 | 1.01e-05 | 2.01e+00
8.84e-03 | 6.24e+03 | 2.41e-06 | 2.07e+00
Table 2: Convergence rates for the truncated fractional diffusion operator (25), $\delta=0.2$ and $h\to 0$ in a continuous Galerkin ansatz space. $h$ | dof | L2 Error | Rates
---|---|---|---
1.41e-01 | 1.44e+02 | 7.58e-04 | 0.00e+00
7.07e-02 | 5.94e+02 | 1.60e-04 | 2.24e+00
3.54e-02 | 2.39e+03 | 4.06e-05 | 1.98e+00
1.77e-02 | 9.59e+03 | 1.01e-05 | 2.01e+00
8.84e-03 | 3.84e+04 | 2.43e-06 | 2.06e+00
Table 3: Convergence rates for the truncated fractional diffusion operator
(25), $\delta=0.2$ and $h\to 0$ in a discontinuous Galerkin ansatz space.

###  3.2 Bond–Based Peridynamics

The translationally invariant and symmetric linear peridynamic kernel is given
by

$\displaystyle\mathbf{C}_{\delta}(\mathbf{x},\mathbf{y}):=c_{\delta}~{}C(\mathbf{x}-\mathbf{y})~{}\mathds{1}_{B_{\delta}(\mathbf{x})}(\mathbf{y}),$

where

$\displaystyle
C(\mathbf{x}-\mathbf{y}):=\frac{(\mathbf{x}-\mathbf{y})(\mathbf{x}-\mathbf{y})^{\top}}{|\mathbf{x}-\mathbf{y}|^{3}}~{}~{}\text{and}~{}~{}c_{\delta}:=\frac{3}{\delta^{3}}.$
(26)

The corresponding linear peridynamic operator then reads as

$\displaystyle-\mathcal{P}_{\delta}\mathbf{u}(\mathbf{x}):=c_{\delta}\int_{B_{\delta}(\mathbf{x})}C(\mathbf{x}-\mathbf{y})(\mathbf{u}(\mathbf{x})-\mathbf{u}(\mathbf{y}))d\mathbf{y}.$
(27)

In [22] the wellposedness is established for problem (21), where
$\mathbf{f}\in L^{2}(\Omega,\mathbb{R}^{d})$ and $\mathbf{g}\in
L^{2}(\Omega_{D},\mathbb{R}^{d})$. For the given constant in (26) it is also
shown there, that the peridynamic operator $-\mathcal{P}_{\delta}$ converges
to the local Navier operator

$\displaystyle-\mathcal{P}_{0}\mathbf{u}(\mathbf{x}):=-\frac{\pi}{4}\Delta\mathbf{u}(\mathbf{x})-\frac{\pi}{2}\nabla\operatorname{div}\mathbf{u}(\mathbf{x})$
(28)

as $\delta\to 0$. A similar convergence result is also obtained for the
corresponding weak solutions. Thus, in the given example we choose the
manufactured polynomial solution
$\mathbf{u}(\mathbf{x}):=(\mathbf{x}_{2}^{2},\mathbf{x}_{1}^{2}\mathbf{x}_{2})$,
set
$\mathbf{f}(\mathbf{x}):=-\mathcal{P}_{0}\mathbf{u}(\mathbf{x})=-\frac{\pi}{2}\left(1+2\mathbf{x}_{1},\mathbf{x}_{2}\right)$
in $\Omega$ and as Dirichlet constraints choose
$\mathbf{g}(\mathbf{x}):=\mathbf{u}(\mathbf{x})$ on $\Omega_{D}$. Similarly to
the diffusion case above, we again obtain that $\mathbf{u}(\mathbf{x})$ is the
solution of (21) due to the correct scaling of the operator [22]. The results
are presented in Table 4 and Table 5 for a continuous and discontinuous
Galerkin ansatz, respectively. In both settings we observe second–order
convergence rate as the mesh size $h\to 0$.

$h$ | dof | L2 Error | Rates
---|---|---|---
1.41e-01 | 3.20e+01 | 7.47e-04 | 0.00e+00
7.07e-02 | 1.62e+02 | 1.82e-04 | 2.04e+00
3.54e-02 | 7.22e+02 | 4.58e-05 | 1.99e+00
1.77e-02 | 3.04e+03 | 1.12e-05 | 2.04e+00
8.84e-03 | 1.25e+04 | 2.63e-06 | 2.09e+00
Table 4: Convergence rates for peridynamic operator (27), $\delta=0.1$ and $h\to 0$ in a continuous Galerkin ansatz space. $h$ | dof | L2 Error | Rates
---|---|---|---
1.41e-01 | 2.88e+02 | 6.64e-04 | 0.00e+00
7.07e-02 | 1.19e+03 | 1.86e-04 | 1.84e+00
3.54e-02 | 4.78e+03 | 4.96e-05 | 1.90e+00
1.77e-02 | 1.92e+04 | 1.25e-05 | 1.99e+00
8.84e-03 | 7.68e+04 | 2.92e-06 | 2.10e+00
Table 5: Convergence rates for peridynamic operator (27), $\delta=0.1$ and
$h\to 0$ in a discontinuous Galerkin ansatz space.

###  3.3 Diffusion with Infinity Ball Truncation

Here, we consider a constant kernel truncated by the infinity normball.
Specifically, we choose the constant to be
$c_{\delta}^{\infty}:=\frac{3}{4\delta^{4}}$, which ensures the convergence to
the local Dirichlet problem for vanishing horizon $\delta\to 0$; see, e.g.,
[31]. The nonlocal operator is then given by

$\displaystyle-\mathcal{L}^{\infty}_{\delta}u(\mathbf{x}):=c_{\delta}^{\infty}\int_{B_{\delta}^{\infty}(\mathbf{x})}(u(\mathbf{x})-u(\mathbf{y}))d\mathbf{y}.$
(29)

The truncation by the infinity normball is implemented without geometric
error, which allows numerical tests of the asymptotic compatibility of the
finite element discretization [30]. For the numerical experiment we choose the
manufactured solution
$u(\mathbf{x})=\sin(4\pi\mathbf{x}_{1})\sin(4\pi\mathbf{x}_{2})$, set
$f(\mathbf{x}):=-\Delta
u(\mathbf{x})=32\pi^{2}\sin(4\pi\mathbf{x}_{1})\sin(4\pi\mathbf{x}_{2})$ in
$\Omega$ and $g(\mathbf{x}):=u(\mathbf{x})$ on $\Omega_{D}$. Note that opposed
to the previous examples, now for a fixed $\delta>0$ the function
$-\mathcal{L}_{\delta}^{\infty}u$ differs from $-\Delta u$, since $u$ is not
chosen to be a polynomial of low degree. Thus the solutions of the nonlocal
and the local Dirichlet problem differ from each other and we can better
observe the asymptotic compatibility of the discretization scheme. We run
tests for a fixed mesh size $h$ and vanishing $\delta$, see Table 6, as well
as for a horizon dependent mesh size $h=\sqrt{2}\delta$ and vanishing
$\delta$, see Table 7. In both cases we observe a second–order convergence
rate as $\delta\to 0$.

$h$ | $\delta$ | L2 Error | Rates
---|---|---|---
8.84e-03 | 2.00e-01 | 2.03e-01 | 0.00e+00
8.84e-03 | 1.00e-01 | 3.98e-02 | 2.35e+00
8.84e-03 | 5.00e-02 | 8.81e-03 | 2.18e+00
8.84e-03 | 2.50e-02 | 2.08e-03 | 2.08e+00
8.84e-03 | 1.25e-02 | 5.17e-04 | 2.01e+00
Table 6: Convergence rates for the infinity normball (29), fixed $h$ and $\delta\to 0$ in a continuous Galerkin ansatz space. $h$ | $\delta$ | L2 Error | Rates
---|---|---|---
1.41e-01 | 2.00e-01 | 2.00e-01 | 0.00e+00
7.07e-02 | 1.00e-01 | 4.01e-02 | 2.32e+00
3.54e-02 | 5.00e-02 | 8.85e-03 | 2.18e+00
1.77e-02 | 2.50e-02 | 2.10e-03 | 2.07e+00
8.84e-03 | 1.25e-02 | 5.17e-04 | 2.03e+00
Table 7: Convergence rates for the infinity normball (29), $\delta=\sqrt{2}h$
and $h,\delta\to 0$ in a continuous Galerkin ansatz space.

###  3.4 Parallel Complexity of the Assembly process

The number of elements in each interaction neighborhood grows quadratically in
2d if the diameter of the elements $h$ is decreased for fixed $\delta$. Thus,
the assembly of the system matrix with retriangulations as described in the
beginning of Section 2.4.1 becomes a costly procedure. Therefore, a matrix
free approach is too expensive and we store the stiffness matrix in a sparse
format. Furthermore, it makes sense to share the work among multiple threads.
The multithreading is implemented using OpenMP [10] and the work is shared by
a partitioning of finite elements as given in Algorithm 2, line 1. Due to the
nonlocality of the operator several threads might need to access identical
entries in the global stiffness matrix at the same time to store their
contribution. OpenMP allows so called critical sections to organize the
manipulation of shared variables. This avoids write conflicts but it would
tremendously slow down the computation. In order to avoid a critical section
during the assembly each thread separately allocates its portion of the global
stiffness matrix. The size of the overlap among the submatrices in the
threads, i.e., the amount of additional memory requirements due to the
parallelization, depends on the nonlocal overlap of the subdomains. We can
therefore reduce the memory requirements by partitioning the domain with metis
[21] instead of a scheduler of OpenMP. The submatrices are finally added
together into a single sparse matrix.

Threads | 1 | 2 | 4 | 8 | 16 | 32 | 64
---|---|---|---|---|---|---|---
Time [s] | 2,153 | 1,169 | 583 | 302 | 166 | 107 | 80
Parallel efficiency | - | 0.92 | 0.92 | 0.89 | 0.81 | 0.63 | 0.42
Table 8: Assembly time for a system with 24,336 degrees of freedom.  Figure
3.1:  Strong scaling for 24,336 (black) and 10,848, 6,120 and 2,736 (gray)
degrees of freedom.

We depict the parallel scaling on a computer with two 2.20GHz Intel Xeon CPUs
with 22 cores per socket and two threads per core. The machine has 756 GB of
RAM. Table 8 and Figure 3.1 (black dots) show the run time of an assembly on a
regular grid on a base domain $\widetilde{\Omega}=[-\delta,0.5+\delta]^{2}$
with mesh size $h=$ 7.1e-03 and $\delta=$ 0.1. The related linear system has
24,336 degrees of freedom and 45,022,167 nonzero entries. The number of
threads is increased by a factor up to 64 while the time drops by a factor of
1/27. The scaling looks perfect for up to 16 parallel threads and the effect
diminishes from then on. Figure 3.1 also shows the scaling for smaller
problems with 10,848, 6,120 and 2,736 degrees of freedom (gray dots) which
show a similar behavior.

##  4 Conclusion

The code nlfem is a flexible tool to set up numerical experiments for
researchers. The documentation also describes the extension by user defined
kernels which allows to consider a large problem class. It can assemble
nonlocal problems in 1d and 3d using the generically implemented barycenter
method [17]. Furthermore, nlfem can be used for a recently developed
substructuring based domain decomposition method for nonlocal operators [9]
which requires reweighting of the kernel function. We therefore hope to bridge
efficiency and flexibility to obtain a convenient Python package which
nourishes the current development in the field of finite element methods for
nonlocal operators and enables easy validations of new theory without the
effort of implementing code from scratch.

## References

  * [1] G. Acosta, F.M. Bersetche, and J.P. Borthagaray.  A short FE implementation for a 2D homogeneous Dirichlet problem of a fractional Laplacian.  Comput. Math. Appl., 74(4):784–816, 2017.
  * [2] M. Ainsworth and C. Glusa.  Towards an efficient finite element method for the integral fractional Laplacian on polygonal domains.  Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan, pages 17–57, 2018.
  * [3] E. Aulisa, G. Capodaglio, A. Chierici, and M. D’Elia.  Efficient quadrature rules for finite element discretizations of nonlocal equations.  Numerical Methods for Partial Differential Equations, 2021.
  * [4] F. Bobaru and W. Hu.  The meaning, selection, and use of the peridynamic horizon and its relation to crack branching in brittle materials.  International Journal of Fracture, 176:215–222, 2012.
  * [5] K. Bogdan, K. Burdzy, and Z.-Q. Chen.  Censored stable processes.  Probability theory and related fields, 127(1):89–152, 2003.
  * [6] B. Boys, T.J. Dodwell, M. Hobbs, and M. Girolami.  Peripy-a high performance opencl peridynamics package.  Computer Methods in Applied Mechanics and Engineering, 386:114085, 2021.
  * [7] D. Brockmann.  Anomalous diffusion and the structure of human transportation networks.  The European Physical Journal Special Topics, 157(1):173–189, 2008\.
  * [8] D. Brockmann and F. Theis.  Money circulation, trackable items, and the emergence of universal human mobility patterns.  IEEE Pervasive Computing, 7(4):28–35, 2008.
  * [9] G. Capodaglio, M. D’Elia, M. Gunzburger, P. Bochev, M. Klar, and C. Vollmann.  A general framework for substructuring-based domain decomposition methods for models having nonlocal interactions.  Numerical Methods for Partial Differential Equations, 2020.
  * [10] L. Dagum and R. Menon.  Openmp: an industry standard api for shared-memory programming.  IEEE computational science and engineering, 5(1):46–55, 1998.
  * [11] M. D’Elia, Q. Du, C. Glusa, X. Tian, and Z. Zhou.  Numerical methods for nonlocal and fractional models.  ACTA Numerica, 29, 2020.
  * [12] M. D’Elia and M. Gunzburger.  The fractional laplacian operator on bounded domains as a special case of the nonlocal diffusion operator.  Comput. Math. Appl., 66:1245–1260, 2013.
  * [13] Q. Du.  Nonlocal Modeling, Analysis, and Computation.  SIAM, 2019.
  * [14] Q. Du, M. Gunzburger, R. B. Lehoucq, and K. Zhou.  Analysis and approximation of nonlocal diffusion problems with volume constraints.  SIAM Review, 54(4):667–696, 2012.
  * [15] Q. Du, M. Gunzburger, R.B. Lehoucq, and K. Zhou.  Analysis and Approximation of Nonlocal Diffusion Problems with Volume Constraints.  SIAM Review, 54(4):667–696, 2012.
  * [16] Q. Du, M. Gunzburger, R.B. Lehoucq, and K. Zhou.  A nonlocal vector calculus, nonlocal volume-constrained problems, and nonlocal balance laws.  Mathematical Models and Methods in Applied Sciences, 23(03):493–540, 2013.
  * [17] M. D’Elia, M. Gunzburger, and C. Vollmann.  A cookbook for approximating Euclidean balls and for quadrature rules in finite element methods for nonlocal problems.  Mathematical Models and Methods in Applied Sciences, 31(08):1505–1567, 2021.
  * [18] G. Da Fies and M. Vianello.  Algebraic cubature on planar lenses and bubbles.  Dolomites Research Notes on Approximation, 5:7–12, 2012.
  * [19] G. Foghem and M. Kassmann.  A general framework for nonlocal Neumann problems.  https://arxiv.org/abs/2204.06793, 2022.
  * [20] G. Gilboa and S. Osher.  Nonlocal operators with applications to image processing.  Multiscale Modeling & Simulation, 7(3):1005–1028, 2009.
  * [21] G. Karypis and V. Kumar.  A fast and high quality multilevel scheme for partitioning irregular graphs.  SIAM Journal on scientific Computing, 20(1):359–392, 1998.
  * [22] T. Mengesha and D. Qiang.  The bond-based peridynamic system with dirichlet-type volume constraint.  Proc. Roy. Soc. Edinburgh Sect. A, 144(1):161–186, 2014.
  * [23] M. Parks, R. Lehoucq, S. Plimpton, and S. Silling.  Implementing peridynamics within a molecular dynamics code.  Computer Physics Communications, 179(11):777–783, 2008.
  * [24] M.L. Parks, P. Seleson, S.J. Plimpton, S.A. Silling, and R.B. Lehoucq.  Peridynamics with lammps: a user guide, v0. 3 beta.  Sandia Report (2011–8253), 3532, 2011.
  * [25] G. Peyré, S. Bougleux, and L. Cohen.  Non-local Regularization of Inverse Problems, pages 57–68.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2008.
  * [26] X. Yin Q. Du, H. Xie.  On some nonlocal models with radially symmetric interaction domains and the effect of domain truncation, 2021.  arXiv:2109.12485.
  * [27] S. Dipierro, X. Ros-Oton, and E. Valdinoci.  Nonlocal problems with Neumann boundary conditions.  Revista Matematica Iberoamericana, 33(2):377–416, 2017.
  * [28] S.A. Sauter and C. Schwab.  Boundary Element Methods, volume 39 of Springer Series in Computational Mathematics.  Springer Berlin Heidelberg, 2011.
  * [29] S.A. Silling.  Reformulation of elasticity theory for discontinuities and long-range forces.  Journal of the Mechanics and Physics of Solids, 48(1):175 – 209, 2000.
  * [30] X. Tian and Q. Du.  Asymptotically compatible schemes and applications to robust discretization of nonlocal models.  SIAM Journal on Numerical Analysis, 52(4):1641–1665, 2014.
  * [31] C. Vollmann.  Nonlocal Models with Truncated Interaction Kernels - Analysis, Finite Element Methods and Shape Optimization.  Doctoral Thesis, Universität Trier, 2019.
  * [32] C.T. Wu, N. Ma, Y. Guo, W. Hu, K. Takada, H. Okada, and K. Saito.  A dynamic ductile failure analysis of shell structures using a nonlocal XFEM method with experimental validation.  Advances in Engineering Software, 123:1–12, 2018.
  * [33] K. Zhou and Q. Du.  Mathematical and Numerical Analysis of Linear Peridynamic Models with Nonlocal Boundary Conditions.  SIAM Journal on Numerical Analysis, 48(5):1759–1780, 2010.