# On immersed boundary kernel functions: a constrained quadratic minimization
perspective

Amneet Pal Singh Bhalla<EMAIL_ADDRESS>Department of Mechanical
Engineering, San Diego State University, San Diego, CA

###### Abstract

In the immersed boundary (IB) approach to fluid-structure interaction
modeling, the coupling between the fluid and structure variables is mediated
using a regularized version of Dirac delta function. In the IB literature, the
regularized delta functions, also referred to IB kernel functions, are either
derived analytically from a set of postulates or computed numerically using
the moving least squares (MLS) approach. Whereas the analytical derivations
typically assume a regular Cartesian grid, the MLS method is a meshless
technique that can be used to generate kernel functions on complex domains and
unstructured meshes. In this note we take a viewpoint that IB kernel
generation, either analytically or via MLS, is a constrained quadratic
minimization problem. The extremization of a constrained quadratic function is
a broader concept than kernel generation and there are well-established
numerical optimization techniques to solve this problem. For example, we show
that the constrained quadratic minimization technique can be used to generate
one-sided (anisotropic) IB kernels and/or to bound their values.

###### keywords:

_immersed boundary method_ , _fluid-structure interaction_ , _meshless
methods_ , _moving least squares method_ , _quadratic programming_

##  1 Introduction

The immersed boundary (IB) method was originally proposed by Peskin to model
blood flow around heart valves [1]. In the IB formulation, the fluid occupying
the entire computational domain is described by Eulerian variables, whereas
the immersed structure that is allowed to freely cut the background Eulerian
grid is described by Lagrangian variables. The presence of the immersed
structure leads to an additional body force term in the momentum equation. The
Lagrangian structure moves with the local fluid velocity and satisfies the no-
slip boundary condition on the fluid-structure interface. The Eulerian-
Lagrangian coupling in the IB formulation is mediated using regularized
integral kernels. Specifically, IB forces defined on the marker points are
_spread_ to the Eulerian grid from the Lagrangian structure, and conversely,
the fluid velocity is _interpolated_ from the Eulerian mesh to the marker
points. The use of an isotropic and regularized IB kernel _diffuses_ the
fluid-structure interface to nearby grid cells and this smearing of the
interface is proportional to the width of the kernel. This also implies that
the (usual) isotropic IB kernels couple the structural degrees of freedom to
fluid variables on both sides of the interface.

In the literature, IB kernels have been derived either analytically or
computed numerically via the moving least squares (MLS) approach. Peskin [2]
used a set of postulates to derive analytical expressions for several IB
kernels such as the four-, five-, and six-point function. These postulates aim
to achieve several desirable features in an IB kernel, such as polynomial
reproducing conditions (also known as the discrete moment conditions), grid
translational invariance property, and compact kernel support. More recently,
Bao et al. [3] used a weaker second moment condition to derive improved
versions of the standard five-, and six-point IB kernels. A weaker second
moment condition allowed the authors to remove the negative tails of the
standard five- and six-point kernels, and to enhance their gird translational
invariance property. A main limitation of the analytical approach is the
assumption of a regular Cartesian grid and the tensor-product form of the
kernel function in dimensions higher than one. To circumvent these
limitations, Vanella and Balaras introduced the MLS approach to generate IB
kernel functions [4]. Although the authors in [4] demonstrated the IB/MLS
approach on a regular Cartesian grid, the technique is easily extendable to
unstructured grid-based IB methods as shown by Krishnan et al. [5] and Saadat
et al. [6]. In this work we take a viewpoint that IB kernel generation is a
constrained quadratic minimization problem and both analytical and MLS
formulations can be described using the language of optimization.

In the prior IB/MLS works [4, 7, 8, 5, 9, 6, 10], a standard formulation of
the moving least squares method was used to transform a cubic spline weight
function into a regularized IB kernel function that satisfies zeroth- and
first-order discrete moment conditions. The standard MLS formulation minimizes
the weighted $L^{2}$-norm of the difference between the interpolant and the
given data values. A lesser known formulation of the moving least squares
problem (in the IB literature) is the Backus-Gilbert theory, which was
developed in the context of geophysics [11]. In the Backus-Gilbert
formulation, polynomial reproducing conditions are used in the MLS problem
definition. However, both formulations ultimately lead to the same IB kernel
function as shown in this work. The connection between the two MLS
formulations was first pointed out by Bos and Šalkauskas in [12] and re-
iterated here in the context of constrained quadratic minimization. We also
show that the constrained quadratic minimization approach can be used to
generate discrete moment conditions satisfying anisotropic IB kernels that
effectively couple structural degrees of freedom to fluid variables on only
one side of the fluid-structure interface. This can help to reduce spurious
feedback forcing and “leakage” that are typically observed in the diffuse-
interface IB models. Furthermore, the constrained quadratic minimization
method can be used to bound the weights of the IB kernels.

In the following sections, we first review the two formulations of the moving
least squares problem and describe them using the language of optimization.
Next, we review the set of postulates used by Peskin to derive the four-point
IB kernel and frame it into the language of constrained quadratic
minimization. Finally, we describe methods to generate one-sided and/or
bounded IB kernels and demonstrate this process using a two-dimensional
example.

##  2 Moving least squares method

###  2.1 The standard formulation

The standard formulation of the moving least squares method considers the
following _approximation_ problem. Given a (column) vector of data values
$\bm{\mathrm{g}}=(g(\bm{\mathrm{x}}_{1}),\ldots,g(\bm{\mathrm{x}}_{N}))^{\intercal}$,
defined on $N$ distinct data sites
$\bm{\mathcal{X}}=\\{\bm{\mathrm{x}}_{1},\ldots,\bm{\mathrm{x}}_{N}\\}$ by a
(smooth) function $g$, find the quasi-interpolant ${\mathcal{P}}g$ to $g$

${\mathcal{P}}g(\bm{\mathrm{x}})=\sum_{j=1}^{m}c_{j}(\bm{\mathrm{x}})p_{j}(\bm{\mathrm{x}})=\bm{\mathrm{p}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}}),$
(1)

at an _evaluation_ point $\bm{\mathrm{x}}$ using the polynomial approximation
space $\mathbb{P}=\\{p_{1},\ldots,p_{m}\\}$. Typically, $m\ll N$ and the
coefficient vector $\bm{\mathrm{c}}(\bm{\mathrm{x}})$ in Eq. (1) is found by
minimizing the weighted $L^{2}$-norm of the error function
$J(\bm{\mathrm{c}})$

$J(\bm{\mathrm{c}}(\bm{\mathrm{x}}))=\sum_{i=1}^{N}W(\bm{\mathrm{x}}_{i},\bm{\mathrm{x}})[\bm{\mathrm{p}}^{\intercal}(\bm{\mathrm{x}}_{i})\bm{\mathrm{c}}(\bm{\mathrm{x}})-g(\bm{\mathrm{x}}_{i})]^{2},$
(2)

with respect to $\bm{\mathrm{c}}$. In the equation above,
$W(\bm{\mathrm{x}})=W(\bm{\mathrm{x}}_{i},\bm{\mathrm{x}})$ is a given weight
function which decreases in magnitude away from the evaluation point and is
taken to be non-negative. The error function $J$ is minimized using the
_essential_ condition of extremum, $\frac{\partial
J}{\partial\bm{\mathrm{c}}}=\bm{\mathrm{0}}$, which leads to the normal
equation of the form

$\sum_{j=1}^{m}c_{j}(\bm{\mathrm{x}})\left<p_{j},p_{k}\right>_{W(\bm{\mathrm{x}})}=\left<g,p_{k}\right>_{W(\bm{\mathrm{x}})},\quad
k=1,\ldots,m.$ (3)

Here, $\left<\cdot,\cdot\right>_{W(\bm{\mathrm{x}})}$ denotes the weighted
inner product of two quantities. If
$\bm{\mathrm{\mathcal{A}}}\in\mathbb{R}^{m\times n}$ denotes the polynomial
matrix with entries
$\bm{\mathrm{\mathcal{A}}}_{ij}=p_{i}(\bm{\mathrm{x}}_{j}),i=1,\ldots,m,j=1,\ldots,N$,
and
$\bm{\mathrm{\mathcal{W}}}(\bm{\mathrm{x}})=\text{diag}\left(\bm{\mathrm{W}}\right)$
is the diagonal matrix containing weights (with
$\bm{\mathrm{W}}=\left[W(\bm{\mathrm{x}}_{1},\bm{\mathrm{x}}),\dots,W(\bm{\mathrm{x}}_{N},\bm{\mathrm{x}})\right]$
as the main diagonal of $\bm{\mathrm{\mathcal{W}}}$), then in matrix notation
Eq. (3) reads as

$\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})=\bm{\mathrm{g}}_{p}(\bm{\mathrm{x}}).$
(4)

Here, $\bm{\mathrm{\mathcal{G}}}$ is the positive-definite Gram matrix
$\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})=\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\mathcal{W}}}\bm{\mathrm{\mathcal{A}}}^{\intercal}\in\mathbb{R}^{m\times
m}$ and $\bm{\mathrm{g}}_{p}$ is the right-hand side vector
$\bm{\mathrm{g}}_{p}(\bm{\mathrm{x}})=\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\mathcal{W}}}\bm{\mathrm{g}}\in\mathbb{R}^{m}$.
Note that the coefficient vector
$\bm{\mathrm{c}}(\bm{\mathrm{x}})=\bm{\mathrm{\mathcal{G}}}^{-1}\bm{\mathrm{g}}_{p}$,
Gram matrix $\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})$, and the right hand
side vector $\bm{\mathrm{g}}_{p}(\bm{\mathrm{x}})$, all depend on the
evaluation point $\bm{\mathrm{x}}$. Consequently, the normal Eq. (4) is solved
anew each time the evaluation point is changed or “moved.”

In the context of the immersed boundary method, the standard MLS formulation
is used to interpolate velocity from Eulerian grid nodes to a Lagrangian
marker point and to spread IB force from a Lagrangian marker point to the
surrounding Eulerian nodes [4]. This is achieved by considering a particular
Lagrangian marker point as the current evaluation point and the surrounding
Eulerian grid nodes (to the Lagrangian point) as the distinct data sites that
provide velocity or force data. The IB/MLS kernel function
$\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\left(\psi_{1}(\bm{\mathrm{x}}),\ldots,\psi_{N}(\bm{\mathrm{x}})\right)^{\intercal}$,
also known as the _generating function_ , is obtained by combining Eqs. (1)
and (3)

$\displaystyle{\mathcal{P}}g(\bm{\mathrm{x}})=\sum_{i=1}^{N}\psi_{i}(\bm{\mathrm{x}})g(\bm{\mathrm{x}}_{i})=\bm{\mathrm{g}}^{\intercal}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})$
$\displaystyle=\bm{\mathrm{p}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})=\bm{\mathrm{g}}^{\intercal}\bm{\mathrm{\mathcal{W}}}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{A}}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{G}}}^{-1}(\bm{\mathrm{x}})\bm{\mathrm{p}}(\bm{\mathrm{x}}).$
(5)

Therefore, from the equation above we get

$\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{\mathcal{W}}}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{A}}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{G}}}^{-1}(\bm{\mathrm{x}})\bm{\mathrm{p}}(\bm{\mathrm{x}}).$
(6)

Comparing the two forms of the quasi-interpolant ${\mathcal{P}}g$ in Eqs. (1)
and (6), it is clear that the latter form employs a kernel function
$\bm{\mathrm{\Psi}}=\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})$ that depends only on
the evaluation point $\bm{\mathrm{x}}$. In contrast, the former form uses a
coefficient vector
$\bm{\mathrm{c}}=\bm{\mathrm{\mathcal{G}}}^{-1}\bm{\mathrm{g}}_{p}$, which
depends on both the evaluation point $\bm{\mathrm{x}}$ and the data values
$\bm{\mathrm{g}}$, i.e.,
$\bm{\mathrm{c}}=\bm{\mathrm{c}}(\bm{\mathrm{x}};\bm{\mathrm{g}})$. For IB
methods, it is therefore more convenient to work with $\bm{\mathrm{\Psi}}$
directly, as it allows the same kernel function to be used in both velocity
interpolation and force spreading operations. It can be shown that using the
same IB kernel function in the two operators conserve energy during
Lagrangian-Eulerian coupling [2].

Although, it is not explicitly enforced in the error minimization of $J$, the
standard MLS method will reproduce the polynomial basis functions
$p_{1},\ldots,p_{m}$, i.e.,

${\mathcal{P}}p_{j}(\bm{\mathrm{x}})=\bm{\mathrm{p}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}};\bm{\mathrm{p}}_{j})=\bm{\mathrm{p}}_{j}^{\intercal}(\bm{\mathrm{x}}_{i})\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=p_{j}(\bm{\mathrm{x}}),\quad\quad\forall
p_{j}\in\mathbb{P}.$ (7)

We shall prove this property of standard MLS formulation in Sec. 2.2.

Next, we show that the standard MLS formulation can be viewed as a constrained
quadratic minimization problem of the form

$\text{Problem
A}=\begin{cases}\min\limits_{\bm{\mathrm{c}}(\bm{\mathrm{x}})\in\mathbb{R}^{m}}&J\;=\;\frac{1}{2}\,\bm{\mathrm{c}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})\\\
\text{subject
to:}&\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})=\bm{\mathrm{g}}_{p}(\bm{\mathrm{x}}).\end{cases}$
(8)

Note that Problem A is a minimization problem because
$\bm{\mathrm{\mathcal{G}}}$ is a (symmetric) positive-definite matrix. The
above constrained optimization problem can be solved by introducing the
Lagrangian function ${\mathcal{L}}$, with the help of a Lagrange multiplier
$\bm{\mathrm{\lambda}}(\bm{\mathrm{x}})\in\mathbb{R}^{m}$

${\mathcal{L}}=\frac{1}{2}\,\bm{\mathrm{c}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})-\bm{\mathrm{\lambda}}^{\intercal}(\bm{\mathrm{x}})[\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})-\bm{\mathrm{g}}_{p}(\bm{\mathrm{x}})].$
(9)

Using the _essential_ conditions of extrema,
$\frac{\partial{\mathcal{L}}}{\partial\bm{\mathrm{c}}}=\bm{\mathrm{0}}$ and
$\frac{\partial{\mathcal{L}}}{\partial\bm{\mathrm{\lambda}}}=\bm{\mathrm{0}}$,
we get the following equalities

$\displaystyle\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})$
$\displaystyle=\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{\lambda}}(\bm{\mathrm{x}}),$
(10)
$\displaystyle\bm{\mathrm{\mathcal{G}}}(\bm{\mathrm{x}})\bm{\mathrm{c}}(\bm{\mathrm{x}})$
$\displaystyle=\bm{\mathrm{g}}_{p}(\bm{\mathrm{x}}),$ (11)

respectively. Solving the above two equations, we obtain
$\bm{\mathrm{c}}(\bm{\mathrm{x}})=\bm{\mathrm{\lambda}}(\bm{\mathrm{x}})=\bm{\mathrm{\mathcal{G}}}^{-1}(\bm{\mathrm{x}})\bm{\mathrm{g}}_{p}(\bm{\mathrm{x}})$.
Hence, the solutions to the standard MLS formulation and Problem A are the
same.

###  2.2 The Backus-Gilbert formulation

The Backus-Gilbert formulation of the moving least squares method seeks the
quasi-interpolant of the type

${\mathcal{P}}g(\bm{\mathrm{x}})=\sum_{i=1}^{N}g(\bm{\mathrm{x}}_{i})\psi_{i}(\bm{\mathrm{x}}),$
(12)

that satisfies the polynomial reproducing conditions

$\sum_{i=1}^{N}p(\bm{\mathrm{x}}_{i})\psi_{i}(\bm{\mathrm{x}})=p(\bm{\mathrm{x}}),\quad\quad\forall
p\in\mathbb{P}.$ (13)

The polynomial reproduction constraints correspond to the _discrete moment
conditions_ for the function $\psi_{i}(\bm{\mathrm{x}})$. Eq. (13) in matrix
form reads as
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$,
in which $\bm{\mathrm{\mathcal{A}}}$ is the same polynomial matrix as defined
in the previous section 2.1. The desired generating function
$\bm{\mathrm{\Psi}}$ in the Backus-Gilbert formulation is obtained by solving
the following constrained quadratic minimization problem:

$\text{Problem
B}=\begin{cases}\min\limits_{\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\in\mathbb{R}^{N}}&J\;=\;\frac{1}{2}\,\bm{\mathrm{\Psi}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{W}}}^{-1}(\bm{\mathrm{x}})\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\\\
\text{subject
to:}&\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}}).\end{cases}$
(14)

The solution to Problem B is obtained by minimizing the Lagrangian function
${\mathcal{L}}$, given by

${\mathcal{L}}=\frac{1}{2}\,\bm{\mathrm{\Psi}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{W}}}^{-1}(\bm{\mathrm{x}})\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})-\bm{\mathrm{\lambda}}^{\intercal}(\bm{\mathrm{x}})[\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})-\bm{\mathrm{p}}(\bm{\mathrm{x}})],$
(15)

in which $\bm{\mathrm{\lambda}}(\bm{\mathrm{x}})\in\mathbb{R}^{m}$ is the
Lagrange multiplier. The Lagrangian ${\mathcal{L}}$ is minimized by using the
essential conditions of extrema,
$\frac{\partial{\mathcal{L}}}{\partial\bm{\mathrm{\Psi}}}=\bm{\mathrm{0}}$ and
$\frac{\partial{\mathcal{L}}}{\partial\bm{\mathrm{\lambda}}}=\bm{\mathrm{0}}$,
which yields the following equalities

$\displaystyle\bm{\mathrm{\mathcal{W}}}^{-1}\bm{\mathrm{\Psi}}-\bm{\mathrm{\mathcal{A}}}^{\intercal}\bm{\mathrm{\lambda}}$
$\displaystyle=\bm{\mathrm{0}},$ (16)
$\displaystyle\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}-\bm{\mathrm{p}}$
$\displaystyle=\bm{\mathrm{0}}.$ (17)

respectively. Solving the above two equations, we obtain

$\displaystyle\bm{\mathrm{\lambda}}$
$\displaystyle=\left(\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\mathcal{W}}}\bm{\mathrm{\mathcal{A}}}^{\intercal}\right)^{-1}\bm{\mathrm{p}}=\bm{\mathrm{\mathcal{G}}}^{-1}\bm{\mathrm{p}},$
(18) $\displaystyle\bm{\mathrm{\Psi}}$
$\displaystyle=\bm{\mathrm{\mathcal{W}}}\bm{\mathrm{\mathcal{A}}}^{\intercal}\bm{\mathrm{\lambda}}=\bm{\mathrm{\mathcal{W}}}\bm{\mathrm{\mathcal{A}}}^{\intercal}\left(\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\mathcal{W}}}\bm{\mathrm{\mathcal{A}}}^{\intercal}\right)^{-1}\bm{\mathrm{p}}=\bm{\mathrm{\mathcal{W}}}\bm{\mathrm{\mathcal{A}}}^{\intercal}\bm{\mathrm{\mathcal{G}}}^{-1}\bm{\mathrm{p}}.$
(19)

Comparing Eqs. (19) and (6), it is clear that both formulations of the moving
least squares method yield the same quasi-interpolant ${\mathcal{P}}g$ to
function $g$. Since the polynomial reproduction constraints are directly
included in the quadratic function minimization used to solve Problem B, the
equivalence of the two formulations implies that the standard moving least
squares method also reproduces the polynomial basis functions.

In the IB literature, constant and linear polynomials are typically employed
in the moving least squares formulation. Defining the polynomials relative to
the evaluation point $\bm{\mathrm{x}}=(\bar{x},\bar{y},\bar{z})$, the set of
basis functions $\mathbb{P}$ for IB problems consists of

$\mathbb{P}=\\{1,x-\bar{x},y-\bar{y},z-\bar{z}\\}.$ (20)

Note that if the weight function $W(\bm{\mathrm{x}}_{i},\bm{\mathrm{x}})$
already satisfies the discrete moment conditions on the support of the IB/MLS
kernel, then the constrained quadratic minimization would result in a trivial
solution $\bm{\mathrm{\Psi}}=\bm{\mathrm{\mathcal{W}}}$. Consequently, the
product $\bm{\mathrm{\mathcal{A}}}^{\intercal}\bm{\mathrm{\lambda}}$ in Eq.
(19) produces a vector whose entries are all one.

##  3 Peskin’s IB kernel

If the IB simulations are performed on a regular Cartesian grid (which is
often the case), then Peskin’s delta functions can be used instead of the
IB/MLS kernels described earlier. Peskin’s delta functions satisfy the
polynomial reproducing conditions, can be expressed compactly in a piecewise
analytical manner, and furthermore, they do not explicitly depend upon a
particular evaluation/marker point. The latter property is due to the fact
that Peskin’s delta functions are constructed in terms of relative distance
between the data site and the evaluation point, i.e,
$\psi=\psi(|\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}|)$ and
$\psi\neq\psi(\bm{\mathrm{x}}_{i},\bm{\mathrm{x}})$.

Taking the example of widely popular four-point delta function $\psi_{4}$ from
the IB literature, we list the four postulates used by Peskin for its
construction below:

$\displaystyle(i)\text{ Zeroth moment:}$
$\displaystyle\quad\quad\sum_{i=1}^{N}\psi_{4}(|\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}|)=1,$
$\displaystyle(ii)\text{ Even-odd:}$
$\displaystyle\quad\quad\sum_{i\;\text{even}}\psi_{4}(|\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}|)=\sum_{i\;\text{odd}}\psi_{4}(|\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}|)=\frac{1}{2},$
$\displaystyle(iii)\text{ First moment:}$
$\displaystyle\quad\quad\sum_{i=1}^{N}(\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}})\psi_{4}(|\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}|)=0,$
$\displaystyle(iv)\text{ Sum-of-squares:}$
$\displaystyle\quad\quad\sum_{i=1}^{N}\psi_{4}^{2}(|\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}|)=\left(\frac{3}{8}\right)^{d}.\quad\quad(d\text{
denotes the spatial dimension})$

Note that the even-odd condition in the second postulate implies the zeroth
moment condition (or the constant reproducing condition), and therefore, the
first equation is not explicitly used in deriving $\psi_{4}$. The first moment
condition in the third postulate is the linear polynomial reproducing
condition, whereas the sum-of-squares condition in the last postulate is a
weaker condition on grid translational invariance.

Similar to the IB/MLS kernel formulation, Peskin’s four-point delta function
can be formulated as a constrained quadratic minimization problem. The
minimization problem reads as

$\text{Problem
C}=\begin{cases}\min\limits_{\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\in\mathbb{R}^{N}}&J\;=\;\frac{1}{2}\,\bm{\mathrm{\Psi}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{W}}}^{-1}(\bm{\mathrm{x}})\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\\\
\text{subject
to:}&(\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}})^{\intercal}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{0}},\\\
&\bm{\mathrm{1}}_{\rm
even}^{\intercal}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\frac{1}{2},\\\
&\bm{\mathrm{1}}_{\rm
odd}^{\intercal}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\frac{1}{2},\\\
&\bm{\mathrm{\Psi}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\Psi}}=\left(\frac{3}{8}\right)^{d},\end{cases}$
(21)

in which $\bm{\mathrm{1}}_{\rm even}$ is a (column) vector containing the
values one and zero at even- and odd-indexed data sites, respectively.
$\bm{\mathrm{1}}_{\rm odd}$ is defined analogously. Note that Problem C has a
nonlinear equality constraint, which would require a nonlinear optimizer, such
as the interior point algorithm [13] or the sequential quadratic programming
technique [14].

##  4 One-sided and bounded IB kernels

Next, we show that the constrained quadratic minimization approach can help to
generate IB kernels that effectively couple structural degrees of freedom to
fluid variables on only one side of the fluid-structure interface. Such
interpolation kernels have been termed “one-sided IB kernels” by Bale et al.
[15], who have recently proposed them. The one-sided IB kernels help to reduce
spurious feedback forcing and internal flows that are typically observed in IB
models that use isotropic kernel functions which couple the structure to fluid
degrees of freedom on both sides of the interface [15]. However, such kernels
have large positive and negative weights, which can lead to numerical
instability over time. To mitigate the numerical instability the authors in
[15] mollified the one-sided IB kernel weights at the expense of the first
moment condition. In contrast, here we show that the numerical values of the
one-sided kernel function can be bounded, e.g., made positive or constrained
within a smaller range using the constrained quadratic minimization approach,
while still satisfying the first moment condition.

(a) Support of an isotropic IB kernel

(b) Support of a one-sided IB kernel

Figure 1: Representative Lagrangian markers (${\color[rgb]{1,0,0}*}$) on a
circular interface embedded in a Cartesian grid.1(a) Full support of a typical
isotropic IB kernel and 1(b) restricted support of one-sided IB kernel. The
Eulerian support of the maker points are shown by circles ($\circ$).

We demonstrate the one-sided IB kernel construction by considering a circular
interface that is embedded in a regular Cartesian grid. The interface is
discretized using Lagrangian marker points, four of which are highlighted by
red markers in Fig. 1. Typical isotropic IB kernels consider both sides of the
interface to define the support region; this is shown by black circles for the
four marker points in Fig. 1(a). For the one-sided kernel, the Eulerian
support region is restricted to only one side of the interface. For example,
in Fig. 1(b) the kernel support is restricted to the exterior of the cylinder,
which is denoted $\Omega^{+}$ region. Analogously, the kernel support can be
defined interior to the cylinder ($\Omega^{-}$ region); here we discuss the
approach in the context of the former case. The one-sided, bounded, and
discrete moment (polynomial reproducing) conditions satisfying IB kernels can
be constructed numerically by solving a constrained quadratic minimization
problem of the type

$\text{Problem
D}=\begin{cases}\min\limits_{\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\in\mathbb{R}^{N}}&J\;=\;\frac{1}{2}\,\bm{\mathrm{\Psi}}^{\intercal}(\bm{\mathrm{x}})\bm{\mathrm{\mathcal{W}}}^{-1}_{\rm
R}(\bm{\mathrm{x}})\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\\\ \text{subject
to:}&\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}},\\\
&\bm{\mathrm{\Psi}}_{\rm
L}\leq\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\leq\bm{\mathrm{\Psi}}_{\rm
U}\end{cases}$ (22)

In the Problem D above, $\bm{\mathrm{\mathcal{W}}}_{\rm R}$ is the restricted
version of the weight function that is defined to be

$\bm{\mathrm{\mathcal{W}}}_{\rm
R}=\begin{cases}W(\bm{\mathrm{x}}_{i},\bm{\mathrm{x}})&\text{ if
}\bm{\mathrm{x}}_{i}\in\Omega^{+},\\\ 0&\text{ if
}\bm{\mathrm{x}}_{i}\in\Omega^{-}.\end{cases}$ (23)

The generating function $\bm{\mathrm{\Psi}}$ in Eq. (22) is bounded on the
lower end by $\bm{\mathrm{\Psi}}_{\rm L}$ values and on the upper end by
$\bm{\mathrm{\Psi}}_{\rm U}$ values. In lack of any particular requirement,
all weights are bounded alike, i.e., $\bm{\mathrm{\Psi}}_{\rm
L}=\alpha\bm{\mathrm{1}}=\bm{\mathrm{\alpha}}$ and $\bm{\mathrm{\Psi}}_{\rm
U}=\beta\bm{\mathrm{1}}=\bm{\mathrm{\beta}}$, in which $\alpha$ and $\beta$
are two scalars and $\bm{\mathrm{1}}$ is a column vector with all entries as
one. Note that we have included an inequality constraint, along with the
equality constraint in Problem D. Without the inequality constraint, Problem D
is same as Problem A or Problem B.

Next, we solve Problem D for the circular interface example by taking constant
and linear polynomial basis functions in $\mathbb{P}$ (as in Eq. (20)). A
linear field of the type

$g(x,y)=10x+5y$ (24)

is defined on the cell centers of a uniform Cartesian mesh that discretizes
the computational domain of extents $\Omega\in[-1,1]^{2}$. The uniform cell
size is taken to be $h=0.075$. As the circular interface partitions the domain
into inner $\Omega^{-}$ and outer $\Omega^{+}$ regions,
$\Omega=\Omega^{-}\cup\Omega^{+}$. The cylinder is centered around the point
$(0,0)$ and has a radius of $R=0.5$. The four marker points on the periphery
of the cylinder make an angle of $40^{\circ}$, $140^{\circ}$ $230^{\circ}$,
and $310^{\circ}$ with the $x$-axis. A tensor-product form of the six-point
spline kernel $\psi_{6}$ is used as a weight function in Problem D. The one-
dimensional version of $\psi_{6}$ function reads as

$\psi_{6}(r)=\begin{cases}\frac{1}{60}(-5\kappa^{5}+90\kappa^{4}-630\kappa^{3}+2130\kappa^{2}-3465\kappa+2193)&0\leq|r|<1,\\\
\frac{1}{120}(5\kappa^{5}-120\kappa^{4}+1140\kappa^{3}-5340\kappa^{2}+12270\kappa-10974)&1\leq|r|<2,\\\
\frac{1}{120}(-\kappa^{5}+30\kappa^{4}-360\kappa^{3}+2160\kappa^{2}-6480\kappa+7776)&2\leq|r|<3,\\\
0,&3\leq|r|.\end{cases}$ (25)

In Eq. (25) above, $r=(x-\bar{x})/h$ and $\kappa=|r|+3$. Note that the weight
function $\psi_{6}$ satisfies constant and linear polynomial reproducing
conditions when the support of the Lagrangian marker is included on both sides
of the interface.

(a) Case 1

(b) Case 2

(c) Case 3

(d) Case 4

Figure 2: Color plot of $\bm{\mathrm{\Psi}}$ for four representative
Lagrangian markers ($\blacktriangle$) placed on a circular interface of radius
$R=0.5$. $\bm{\mathrm{\Psi}}$ is obtained using the constrained quadratic
minimization formulation of Problem D. The interface is embedded in a
computational domain of extents $\Omega\in[-1,1]^{2}$. Interpolation weights
when 2(a) Eulerian support on both sides of the interface and only equality
constraints
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
are used; 2(b) Eulerian support restricted to $\Omega^{+}$ region and only
equality constraints
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
are used; 2(c) Eulerian support restricted to $\Omega^{+}$ region and both
equality
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
and inequality constraints
$\bm{\mathrm{-0.07}}\leq\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\leq\bm{\mathrm{0.5}}$
are used; and 2(d) Eulerian support restricted to $\Omega^{+}$ region and both
equality
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
and inequality constraints
$\bm{\mathrm{0}}\leq\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\leq\bm{\mathrm{0.75}}$
are used.

The linear field $g$ is interpolated onto the four marker points using

  * 1.

Case 1: Eulerian support on both sides of the interface and only equality
constraints
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
in Problem D;

  * 2.

Case 2: Eulerian support restricted to $\Omega^{+}$ region and only equality
constraints
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
in Problem D;

  * 3.

Case 3: Eulerian support restricted to $\Omega^{+}$ region and both equality
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
and inequality constraints
$\bm{\mathrm{-0.07}}\leq\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\leq\bm{\mathrm{0.5}}$
in Problem D; and

  * 4.

Case 4: Eulerian support restricted to $\Omega^{+}$ region and both equality
$\bm{\mathrm{\mathcal{A}}}\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})=\bm{\mathrm{p}}(\bm{\mathrm{x}})$
and inequality constraints
$\bm{\mathrm{0}}\leq\bm{\mathrm{\Psi}}(\bm{\mathrm{x}})\leq\bm{\mathrm{0.75}}$
in Problem D.

The three cases listed above are solved using the quadprog routine in MATLAB
using the default solver options. The quadprog routine requires defining the
cost function $J$ in terms of a matrix ($\bm{\mathrm{\mathcal{W}}}^{-1}$ in
this case) and lower and upper bound vectors ($\bm{\mathrm{\Psi}}_{\rm L}$ and
$\bm{\mathrm{\Psi}}_{\rm R}$, respectively in this case) for the variable to
be minimized ($\bm{\mathrm{\Psi}}$ in this case).

Table 1 lists the relative error in the interpolated values of $g$ and Fig. 2
shows the color plot of the interpolation weights around the four marker
points. Fig. 2(a), which plots the weights for Case 1 shows that Problem D
generates the usual isotropic IB kernel because the support of the marker
points is defined on both sides of the interface. In fact for Case 1, Problem
D did not modify the underlying weight function $\bm{\mathrm{\Psi}}_{6}$
because the weight function already satisfies the polynomial reproducing
conditions, as mentioned earlier. This can be verified by computing the
$L^{2}$\- and $L^{\infty}$-norms of the difference
$||\bm{\mathrm{\Psi}}-\bm{\mathrm{\Psi}}_{6}||$, which come out to be
$2.5829\times 10^{-13}$ and $1.0819\times 10^{-13}$, respectively.

Interpolation weights for Case 2 are shown in Fig. 2(b). In this case we
observe large positive weights near the marker points, which are followed by
large negative values; $\bm{\mathrm{\Psi}}$ ranges from -0.3627 to 0.9178 in
Case 2. This kind of weight distribution (having large positive and negative
values) was shown to have a destabilizing effect when used within a direct-
forcing IB method [15]. To make use of the one-sided IB kernel, the authors in
[15] mollified the weights through a renormalization procedure. However, the
mollified weights satisfy only the zeroth moment condition, i.e, they are only
first-order accurate.

Next, we show that Problem D can be used to systematically mollify the non-
smooth weights obtained in Case 2 by prescribing the upper and lower bounds,
without having to forfeit the first moment condition. The resulting smooth
weights for Case 3 are shown in Fig. 2(c), wherein we take $\alpha=-0.07$ and
$\beta=0.5$, which is a subset of the range observed for $\bm{\mathrm{\Psi}}$
in Case 2. To quantitatively access the accuracy of the IB kernels, we
tabulate the interpolation error $\mathcal{E}$ at the four marker locations in
Table 1. From the table, it can be inferred that close to machine precision
error norms are obtained for Case 1 and 2, which indicates that the linear
field is interpolated exactly onto the Lagrangian marker points. This is not a
surprise, because the linear polynomial reproducing constraint is imposed in
the problem formulation itself, but it demonstrates the procedure to obtain
discrete moment conditions satisfying one-sided kernels using quadratic
programming, which is different compared to the moving-least squares technique
of Problem A or B. For Case 3, $\mathcal{E}$ is in the range of
$10^{-10}-10^{-11}$. We remark that close to machine precision error norms can
be achieved in this case by further decreasing (increasing) the lower (upper)
bound of $\bm{\mathrm{\Psi}}$. However, this will not lead to any practical
difference in the results of an IB simulation and loosing some digits of
interpolation accuracy should be preferred over accurate but non-smooth
weights, such as those shown in Fig. 2(b).

Lastly, we consider Case 4, in which our objective is to completely eliminate
the negative weights of Case 2. To achieve this, we take $\alpha=0$ and
$\beta=0.75$. Table 1 shows the error norms for this case. Here, we observe
that although the error norms are small (in the range of $10^{-7}-10^{-10}$),
they are not close to machine precision, indicating that the linear field is
not interpolated exactly. This is because the quadprog routine is unable to
find weights that can satisfy the equality constraints exactly, while
simultaneously respecting the imposed bounds; changing the solver settings did
not affect the final convergence of quadprog (data not shown). On the other
hand, if we observe the min/max limits of the color bar in Fig. 2(d),
$\bm{\mathrm{\Psi}}$ for Case 4 stays within the imposed limits. Note that
although the one-sided weights of Case 4 are positive, they are not smooth and
are highly localized near the interface. This type of weight distribution can
also lead to numerical instability for the direct-forcing IB method [16].
Apparently, the weight-shifting technique of Bale et al. to completely
eliminate the negative weights at the expense of first moment condition is
more suitable for the direct-forcing IB approach than the inequality
constraint employed in Case 4. Nevertheless, if we allow small magnitude of
negative weights in the constrained quadratic formulation, we can achieve
smoothness and as well as the desired accuracy in the one-sided IB kernels
(Case 3).

Table 1: Relative error $\mathcal{E}=\frac{\left|\mathcal{P}g-g\right|}{|g|}$ for the three cases of Problem D. For Case 3, $\alpha=-0.07$ and $\beta=0.5$ and for Case 4, $\alpha=0$ and $\beta=0.75$ Case | Relative error $(\mathcal{E})$ at marker location
---|---
40∘ | 140∘ | 230∘ | 310∘
1 | $1.6335\times 10^{-16}$ | $1.9975\times 10^{-16}$ | $6.9267\times 10^{-16}$ | $1.1967\times 10^{-15}$
gray!10 2 |  gray!10 $7.6776\times 10^{-15}$ |  gray!10 $2.3970\times 10^{-14}$ |  gray!10 $3.117\times 10^{-15}$ |  gray!10 $5.2997\times 10^{-15}$
3 | $6.1497\times 10^{-11}$ | $1.3954\times 10^{-10}$ | $2.2942\times 10^{-11}$ | $4.8383\times 10^{-10}$
gray!10 4 |  gray!10 $8.2015\times 10^{-11}$ |  gray!10 $2.4047\times 10^{-9}$ |  gray!10 $1.2423\times 10^{-7}$ |  gray!10 $3.0356\times 10^{-9}$

##  5 Summary and Conclusions

In this note we took a viewpoint that IB kernel generation, either
analytically or via MLS, is a constrained quadratic minimization problem. The
extremization of a constrained quadratic function is a broader concept than
kernel generation and there are well-established numerical optimization
techniques to solve this problem. We also reviewed the two formulations of the
moving least squares problem and described them using the language of
optimization. Further, the set of postulates used by Peskin to derive the
four-point IB kernel was also described using this language. Finally, we
described methods to generate one-sided and/or bounded IB kernels and
demonstrated the process using a two-dimensional example. Four constrained
quadratic minimization problems were described, namely, the Problem A, B, C,
and D. All but Problem C involved linear constraints and can be solved using
efficient optimization techniques such as the quadratic programming technique.
For Problem C, a nonlinear optimizer, such as the interior point algorithm or
the sequential quadratic programming technique is required.

## Acknowledgements

A.P.S.B acknowledges helpful discussion with Rahul Bale.

## Bibliography

## References

  * [1] C. S. Peskin, Numerical analysis of blood flow in the heart, Journal of computational physics 25 (3) (1977) 220–252.
  * [2] C. S. Peskin, The immersed boundary method, Acta Numer 11 (2002) 479–517.
  * [3] Y. Bao, J. Kaye, C. S. Peskin, A gaussian-like immersed-boundary kernel with three continuous derivatives and improved translational invariance, Journal of Computational Physics 316 (2016) 139–144.
  * [4] M. Vanella, E. Balaras, Short note: A moving-least-squares reconstruction for embedded-boundary formulations, Journal of Computational Physics 228 (18) (2009) 6617–6628.
  * [5] S. Krishnan, E. S. Shaqfeh, G. Iaccarino, Fully resolved viscoelastic particulate simulations using unstructured grids, Journal of Computational Physics 338 (2017) 313–338.
  * [6] A. Saadat, C. J. Guido, G. Iaccarino, E. S. Shaqfeh, Immersed-finite-element method for deformable particle suspensions in viscous and viscoelastic media, Physical Review E 98 (6) (2018) 063316.
  * [7] D. Li, A. Wei, K. Luo, J. Fan, An improved moving-least-squares reconstruction for immersed boundary method, International Journal for Numerical Methods in Engineering 104 (8) (2015) 789–804.
  * [8] M. D. de Tullio, G. Pascazio, A moving-least-squares immersed boundary method for simulating the fluid–structure interaction of elastic bodies with arbitrary thickness, Journal of Computational Physics 325 (2016) 201–225.
  * [9] D.-V. Le, B.-C. Khoo, A moving-least-square immersed boundary method for rigid and deformable boundaries in viscous flow, Communications in Computational Physics 22 (4) (2017) 913–934.
  * [10] M. Haji Mohammadi, F. Sotiropoulos, J. Brinkerhoff, Moving least squares reconstruction for sharp interface immersed boundary methods, International Journal for Numerical Methods in Fluids 90 (2) (2019) 57–80.
  * [11] G. Backus, F. Gilbert, The resolving power of gross earth data, Geophysical Journal International 16 (2) (1968) 169–205.
  * [12] L. Bos, K. Salkauskas, Moving least-squares are backus-gilbert optimal, Journal of Approximation Theory 59 (3) (1989) 267–275.
  * [13] A. Wächter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming, Mathematical programming 106 (1) (2006) 25–57.
  * [14] J. Nocedal, S. J. Wright, Sequential quadratic programming, Numerical optimization (2006) 529–562.
  * [15] R. Bale, A. P. S. Bhalla, B. E. Griffith, M. Tsubokura, A one-sided direct forcing immersed boundary method using moving least squares, Journal of Computational Physics 440 (2021) 110359.
  * [16] M. Uhlmann, An immersed boundary method with direct forcing for the simulation of particulate flows, Journal of Computational Physics 209 (2) (2005) 448–476.