STATUS:  Prepared

(1.25in,2in) Numerical simulation of an extensible capsule using regularized
Stokes kernels and overset finite differences

Dhwanit Agarwal
_Oden Institute, University of Texas at Austin
Austin, TX, 78712_
<EMAIL_ADDRESS>
George Biros
_Oden Institute, University of Texas at Austin
Austin, TX, 78712_
<EMAIL_ADDRESS>
(1.25in,7in)

###### Abstract

In this paper, we present a novel numerical scheme for simulating deformable
and extensible capsules suspended in a Stokesian fluid. The main feature of
our scheme is a partition-of-unity (POU) based representation of the surface
that enables asymptotically faster computations compared to spherical-
harmonics based representations. We use a boundary integral equation
formulation to represent and discretize hydrodynamic interactions. The
boundary integrals are weakly singular. We use the quadrature scheme based on
the regularized Stokes kernels (given in [32]). We also use partition-of unity
based finite differences that are required for the computational of
interfacial forces. Given an $N$-point surface discretization, our numerical
scheme has fourth-order accuracy and $\mathcal{O}(N)$ asymptotic complexity,
which is an improvement over the $\mathcal{O}(N^{2}\log N)$ complexity of a
spherical harmonics based spectral scheme that uses product-rule quadratures
[34]. We use GPU acceleration and demonstrate the ability of our code to
simulate the complex shapes with high resolution. We study capsules that
resist shear and tension and their dynamics in shear and Poiseuille flows. We
demonstrate the convergence of the scheme and compare with the state of the
art.

###### Contents

  1. 1 Introduction
  2. 2 Problem formulation
    1. 2.1 Formulation
    2. 2.2 Boundary integral formulation
  3. 3 Numerical algorithms
    1. 3.1 Surface parameterization
    2. 3.2 Surface discretization
    3. 3.3 Smooth surface integrals
    4. 3.4 Singular integration
    5. 3.5 Surface derivatives
    6. 3.6 Time stepping and the overall algorithm
  4. 4 Results
    1. 4.1 Integration results
    2. 4.2 Derivative results
    3. 4.3 Accuracy and convergence of the full numerical scheme
    4. 4.4 Relaxation of capsule
    5. 4.5 Capsule in shear flow
    6. 4.6 Capsule in Poiseuille flow
    7. 4.7 Timing results
  5. 5 Fast multipole method based acceleration
  6. 6 Conclusions
  7. A Appendix
    1. A.1 Transition maps
    2. A.2 Surface derivative formulas
    3. A.3 The choice of parameter $r_{0}$ for the partition of unity
    4. A.4 Effect of the blending process in calculating derivatives
    5. A.5 Surface derivative and singular integration errors for shear flow and Poiseuille flow terminal shapes
    6. A.6 Numerical errors for different reduced volume $\nu$
    7. A.7 Sensitivity of the singular quadrature scheme to the regularization parameter $\delta$
    8. A.8 Wall clock times for our GPU accelerated code _vs_ the spherical harmonics CPU code

##  1 Introduction

Capsules suspended in a Stokesian fluid describe complex biological flows and
microfluidic devices [35, 12, 20]. In particular, we’re interested in modeling
red blood cell (RBC) suspensions. Several previous works [13, 19, 8, 23, 42]
have modeled red blood cells as elastic membranes filled with a Newtonian
fluid and suspended in a Newtonian fluid, also referred to as _“capsules”_. In
this work, we present a numerical scheme for simulating a single deformable
three-dimensional capsule suspended in Stokes flow where its membrane resists
shear and tension [30, 29] (see Figure 1). The scheme can be extended using
existing techniques to simulating an arbitrary number of capsules [34].

_Related work:_ Stokesian flows with deformable capsules involve moving
interfaces, fluid-structure interaction, large deformations, near and long-
range stiff interactions, and often require high-order surface derivatives.
Efficient methods for simulating Stokesian capsule suspensions include
immersed boundary/interface methods [3], lattice-Boltzmann methods [22], and
boundary integral equations methods [27, 42, 34, 11]. Each of these methods
trades-off aspects of accuracy, simplicity of implementation, the ability to
simulate complex physics, and the numerical efficiency. Without advocating a
particular approach, we focus on integral equation formulations and discuss
the related literature in somehow greater detail.

A nice feature of boundary integral formulations for Stokes flows is that they
require only the discretization of the capsule surface. One difficulty is that
they require specialized singular quadrature schemes to compute layer
potentials Spectrally accurate boundary integral methods [34, 42] have used
spherical harmonics representations for the surface and surface fields,
_e.g.,_ elastic forces, velocity, and shape, to simulate the time evolution of
the surface. This surface representation is spectrally accurate in the degree
$p$ of the spherical harmonics that translates to $O(p^{2})$ discretization
points. For spectrally accurate singular quadrature, the authors in [42] used
the Bruno-Kunyansky quadrature [7, 41] that uses a floating partition of unity
to resolve weakly singular points around each global quadrature point. This
scheme has $O(p^{4}\log{p})$ work complexity and can be reduced to $O(p^{3})$
using an FFT-based scheme that embeds the surface [7] in a volumetric
Cartesian grid. The authors in [34] adapted the Graham-Sloan quadrature [17]
to the Stokes kernel. Given $\mathcal{O}(p^{2})$ global points to represent
the capsule, the scheme uses a product quadrature rule that requires
$\mathcal{O}(p^{2})$ spherical harmonics rotations. For small $p$, the fastest
implementation uses dense linear algebra at a $\mathcal{O}(p^{6})$ cost; this
becomes prohibitively expensive for large scheme. Fast rotation schemes based
on Legendre transforms and FFTs can lower the cost to $O(p^{4}\log p)$. The
scheme has excellent accuracy, and even for small values of $p$ for nearly-
spherical shapes it outperforms the Bruno-Kunyansky quadrature. This is due to
Bruno-Kunyansky’s highly localized floating-partition of unity that increases
the resolution requirements. Despite the spectral convergence of these methods
in the degree $p$ of spherical harmonics used, the spherical harmonics
representation has difficulties in resolving shapes with high curvature, like
the shapes shown in Figure 13(c) and Figure 16(c). Indeed, increasing the
degree $p$ of the spherical harmonics representation is too expensive [1, 2]
because they are not amenable to acceleration via fast multipole methods
because of the product quadrature. An alternative is to use more classical
surface representation that avoids global polynomials. Several studies [6, 5,
43, 31, 11] have used triangulation based surface representation to provide
more control over local resolution. These schemes are effective but they only
provide second-order accurate representation for the surfaces. Using a lower-
order representation makes the computation of higher order derivatives more
difficult and complicate the quadrature rules, although recent advances in the
latter may circumvent this difficulty [36] and [18]. We remark, that we prefer
schemes that no expensive geometry-dependent rule precomputation as the
capsule geometry evolves in time.

_Our Contributions:_ We present a new high-order numerical scheme for
simulating the dynamics of three-dimensional extensible capsules. The new
scheme allows non-uniform discretization and fast evaluations of integral
operators. Our singular quadrature is based on the regularized Stokes kernels
in [32]. We use multiple overlapping patches to parameterize the capsule
surface assuming a spherical topology surface at all times. We use a fourth-
order finite difference scheme using an overset grid based discretization of
the patches to calculate the interfacial elastic forces. For $N$ number of
discretization points of the surface, our scheme has $O(N^{2})$ work
complexity, which is similar to the lower order boundary element schemes [11]
and is fourth order accurate. Hence, our scheme offers a high order accuracy
and independent control over local resolution. Our singular quadrature scheme
is not a product quadrature, which allows acceleration via fast multipole
methods (FMMs) [33]. With FMM acceleration, our numerical scheme becomes an
$O(N)$ scheme.

To summarize, our scheme comprises the following main components:

  * •

We describe an overlapping patch based parameterization for capsules that are
diffeomorphic to the unit sphere.

  * •

We provide a finite difference scheme based on the overset grid based
discretization of these patches to calculate the surface derivatives.

  * •

We evaluate a high order singular quadrature scheme based on the regularized
Stokes kernels given in [32].

  * •

We use GPU acceleration for the evaluation of singular quadrature to do high
resolution simulations.

_Limitations:_ We assume spherical-topology of the capsule surface which may
not be the case in certain biophysical flows. In this paper, we do not
consider the bending elastic forces [28, 11, 34] but the can be extended as
is. We do not consider the viscosity contrast, _i.e.,_ the fluid inside and
outside the capsule have the same viscosity. Our algorithm can be extended to
include differences in the interior and exterior fluid viscosity by including
a double layer Stokes potential term in the problem formulation as described
in detail in [28]. The regularized Stokes double layer potential is provided
in [32]. Our scheme requires us to choose several parameters that can affect
the accuracy of the scheme, for example, the configuration of patches, the
extent of overlap of the patches, the partition of unity functions, and a
regularization parameter for the singular quadrature. That is, instead of a
regular grid, we eventually want to switch to quadtree-based representation.
We also do not provide curvature adaptive surface parameterization. Our scheme
does not provide an arbitrary order of accuracy. We remark however, that a
recently developed scheme [39] can be used instead of [32] and provide local
quadrature with high accuracy.

_Outline of the paper:_ In the next section, we state the differential and the
boundary integral formulation of the dynamics of an extensible capsule
suspended in Stokes flow. In Section 3, we describe in detail the surface
parameterization and the numerical schemes (differentiation, integration and
time stepping) we use to simulate the capsule dynamics. In Section 4, we
report the numerical results demonstrating the convergence and accuracy of the
numerical schemes along with the results of capsule dynamics in shear and
Poiseuille flow. In Section 5, we discuss the FMM based acceleration of our
scheme and show results for the speedup obtained via acceleration using a
single level FMM. In Section 6, we provide a brief summary of the paper along
with the conclusions and the future directions that can be taken to improve
upon this work.

Figure 1: A representation of the problem setup. The grey filled region is the
interior of the capsule with membrane $\gamma$. Exterior of the capsule is
filled with a Newtonian fluid and the capsule is suspended freely in it.
$\bm{u}_{\infty}$ is the imposed background fluid velocity.

##  2 Problem formulation

In this section, we formally describe the mathematical formulation of the
problem and introduce the notation. We discuss the differential formulation in
Section 2.1 and then specify the boundary integral formulation of the problem
in Section 2.2. A representation of the setup is shown in Figure 1. Our
discussion follows the work in [27, 34].

###  2.1 Formulation

Let $\gamma$ be the membrane of an extensible capsule filled with a viscous
Newtonian fluid of viscosity $\mu_{i}$ and suspended in another viscous
Newtonian fluid of viscosity $\mu_{e}$. In this work, we assume that the
interior and exterior fluids have same viscosity, _i.e.,_
$\mu_{i}=\mu_{e}=\mu$. The microscopic length scale of the problem implies
that the dynamics of the problem can be described by the Stokes equation. The
boundary value problem is given by [34]:

$\displaystyle-\mu\Delta\bm{u}(\bm{x})+\nabla p(\bm{x})$ $\displaystyle=0\
\quad\forall\bm{x}\in\mathbb{R}^{3}\backslash\gamma,$ (2.1)
$\displaystyle\nabla\cdot\bm{u}(\bm{x})$
$\displaystyle=0\quad\forall\bm{x}\in\mathbb{R}^{3}\backslash\gamma,$ (2.2)
$\displaystyle[[-p\bm{n}+\mu(\nabla\bm{u}+\nabla\bm{u}^{T})\bm{n}]]$
$\displaystyle=\bm{f}\quad\textit{ on }\gamma,$ (2.3)
$\displaystyle\bm{u}(\bm{x})$
$\displaystyle\longrightarrow\bm{u}_{\infty}\quad\textit{ for
}\bm{x}\longrightarrow\infty,$ (2.4)
$\displaystyle\frac{\partial\bm{x}}{\partial t}$
$\displaystyle=\bm{u}\quad\textit{ on }\gamma.$ (2.5)

Here $\mu$ is the viscosity of the exterior and interior fluid, $\bm{u}$ is
the velocity of the fluid, $p$ is the pressure, [[$q$]] represents the jump in
a quantity $q$ across the capsule membrane $\gamma$, $\bm{n}$ is the unit
normal to the capsule membrane, $\bm{f}$ is the total force exerted by the
capsule membrane onto the fluid, $\bm{u}_{\infty}$ is the imposed background
velocity far away from the capsule. Equation 2.1 is the Stokes equation
representing the conservation of momentum, Equation 2.2 is the conservation of
mass, Equation 2.3 is the balance of force on the interface between fluid and
membrane, Equation 2.4 sets the far field velocity to be the background
velocity and Equation 2.5 enforces the no-slip condition on the capsule
membrane.

_Interfacial force:_ Now we discuss in detail the interfacial force $\bm{f}$
exerted by capsule membrane onto the surrounding fluid due to membrane
elasticity. We only consider the elastic forces due to in-plane shear
deformations of the capsule. Previous works have also included bending forces
[42, 11, 38] but we leave it for future work and focus mainly on the numerical
schemes in this paper. The in-plane shear force $\bm{f}_{s}$ is equal to the
surface divergence of the symmetric part of the in-plane shear stress tensor
denoted by $\Lambda$ [29, 42]. Thus, we write,

$\displaystyle\bm{f}=\bm{f}_{s}=\nabla_{\gamma}\cdot\Lambda\textit{. }$ (2.6)

_In-plane shear stress tensor $\Lambda$:_ Our discussion follows the work
based on the Skalak model in [27, 29]. The in-plane shear stress tensor
$\Lambda$ is a function of the surface deformation gradient of the capsule
surface relative to the stress-free reference configuration of the membrane.
Let $\gamma$ be the surface in the current configuration and $\gamma_{r}$ be
the reference configuration of the membrane. If $\bm{x}_{r}\in\gamma_{r}$ is a
point that maps to $\bm{x}\in\gamma$ in the current configuration, let
$\xi(\bm{x}_{r})=\bm{x}$ be the bijective map between the configurations. The
deformation gradient $\bm{F}$ is defined as
$\bm{F}=\frac{\partial\xi}{\partial\bm{x}_{r}}.$ If $\bm{n}_{r}$ is the normal
to the reference configuration at $\bm{x}_{r}\in\gamma_{r}$ and $\bm{n}$ is
the normal to the current configuration at $\bm{x}\in\gamma$, then the
relative surface deformation gradient $\bm{F}_{S}$ is defined as
$\bm{F}_{S}=(\bm{I}-\bm{n}\bm{n}^{T})\bm{F}(\bm{I}-\bm{n}_{r}\bm{n}_{r}^{T})$,
where $\bm{I}$ is the identity tensor. It follows from the above equation that
$\bm{F}_{S}$ maps the surface tangents $\bm{a}_{1r}$ and $\bm{a}_{2r}$ in the
reference configuration $\gamma_{r}$ to the tangents $\bm{a}_{1}$ and
$\bm{a}_{2}$ in the current configuration $\gamma$. Also, it maps the
reference normal $\bm{n}_{r}$ to $0$. Thus, the following set of equations at
every point $\bm{x}_{r}$ uniquely determine $\bm{F}_{S}$ at that point on the
reference configuration: $\bm{F}_{S}\bm{a}_{1r}=\bm{a}_{1},\textit{
}\bm{F}_{S}\bm{a}_{2r}=\bm{a}_{2},\textit{ }\bm{F}_{S}\bm{n}_{r}=0.$ The left
Cauchy-Green deformation tensor $\bm{V}^{2}$ and the surface projection tensor
$\bm{P}$ at point $\bm{x}$ in the current configuration are then defined as

$\displaystyle\bm{V}^{2}(\bm{x})=\bm{F}_{s}(\xi^{-1}(\bm{x}))(\bm{F}_{s}(\xi^{-1}(\bm{x})))^{T},\textit{
}\bm{P}(\bm{x})=\bm{I}-\bm{n}\bm{n}^{T}.$ (2.7)

Following the work in [30, 29], if $\lambda_{1}^{2},\lambda_{2}^{2}$ are the
two non-zero eigenvalues of $\bm{V}^{2}$, we define two scalar invariants
$I_{1}=\lambda_{1}^{2}+\lambda_{2}^{2}-2\textit{ and
}I_{2}=\lambda_{1}^{2}\lambda_{2}^{2}-1$. We now define the shear stress
tensor $\Lambda$ at every point $\bm{x}$ on the current configuration as

$\displaystyle\Lambda(\bm{x})=\frac{E_{s}}{2J_{s}}(I_{1}+1)\bm{V}^{2}(\bm{x})+\frac{J_{s}}{2}(E_{D}I_{2}-E_{s})\bm{P}(\bm{x}),$
(2.8)

where $J_{s}=\lambda_{1}\lambda_{2}$, $E_{s}$ is the elastic shear modulus of
the membrane and $E_{D}$ is the dilatation modulus. High values of the shear
modulus $E_{s}$ represent higher membrane shear resistance. The dilatation
modulus $E_{D}$ controls the local membrane extensibility.

###  2.2 Boundary integral formulation

Following the work in [29, 27], Equations 2.1, 2.2, 2.3, 2.4 and 2.5 can be
formulated into integral equations on the capsule membrane $\gamma$. This
formulation can be written as:

$\displaystyle\bm{f}(\bm{x})=\nabla_{\gamma}\cdot\Lambda(\bm{x})\quad\forall\bm{x}\in\gamma,$
(2.9)
$\displaystyle\bm{u}(\bm{x})=\bm{u}_{\infty}(\bm{x})+\mathcal{S}_{\gamma}[\bm{f}](\bm{x})\quad\forall\bm{x}\in\gamma,$
(2.10) $\displaystyle\frac{\partial\bm{x}}{\partial
t}=\bm{u}(\bm{x})\quad\forall\bm{x}\in\gamma.$ (2.11)

Here $\tau(\bm{x})$ is computed using Equation 2.8.
$\mathcal{S}_{\gamma}[\bm{f}]$ represents the single layer potential of layer
density $\bm{f}$ over the capsule membrane $\gamma$ defined as follows:

$\displaystyle\mathcal{S}_{\gamma}[\bm{f}](\bm{x})=\int_{\gamma}\mathcal{G}(\bm{x},\bm{y})\bm{f}(\bm{y})d\gamma,$
(2.12)

where $\mathcal{G}$ is the Stokes kernel given by
$\mathcal{G}(\bm{x},\bm{y})=\frac{1}{8\pi\mu}\left(\frac{I}{||\bm{r}||}+\frac{\bm{r}\otimes\bm{r}}{||\bm{r}||^{3}}\right)$
with $\bm{r}=\bm{x}-\bm{y}$. Given the initial position of the capsule and its
stress-free reference configuration, we will use Equations 2.9, 2.10 and 2.11
to simulate the time evolution of capsule under the imposed background flow
$\bm{u}_{\infty}$.

##  3 Numerical algorithms

In this section, we describe the numerical algorithms we use to simulate the
dynamics of the capsule. We discretize Equations 2.9, 2.10 and 2.11 discussed
in previous section. In Section 3.1, we discuss our parameterization of
capsule surface (assuming it is diffeomorphic to the unit sphere) using
multiple patches. In Section 3.2, we discuss the discretization of the surface
based on our parameterization of the surface. In Section 3.3, we discuss the
numerical scheme for the surface integration of a smooth function. In Section
3.4, we discuss the numerical scheme for the singular integration to compute
Stokes single layer potential. In Section 3.5, we discuss the numerical scheme
for calculating the surface derivatives for the computation of the elastic
force. Finally in Section 3.6, we discuss the time stepping we use to simulate
the capsule dynamics and provide a summary of the overall algorithm.

###  3.1 Surface parameterization

We assume the capsule membrane $\gamma$ to be smooth and diffeomorphic to the
unit sphere $\mathbb{S}^{2}$ embedded in $\mathbb{R}^{3}$. Thus, we have a
smooth bijective map $\phi:\mathbb{S}^{2}\longrightarrow\gamma$ with the
smooth inverse $\phi^{-1}$. We define $n_{p}$ number of overlapping patches
$\\{\mathcal{P}_{i}^{0}\\}_{i=1}^{i=n_{p}}$, where
$\mathcal{P}_{i}^{0}\subset\mathbb{S}^{2}$ form an open cover of the unit
sphere $\mathbb{S}^{2}$, _i.e.,_
$\bigcup_{i=1}^{n_{p}}\mathcal{P}_{i}^{0}=\mathbb{S}^{2}$. Additionally, we
assume each patch $\mathcal{P}_{i}^{0}$ is diffeomorphic to an open set
$\mathcal{U}_{i}$ ($\subset\mathbb{R}^{2}$) via a coordinate chart
$\eta_{i}^{0}$, _i.e.,_
$\eta_{i}^{0}:\mathcal{U}_{i}\longrightarrow\mathcal{P}_{i}^{0}$ is a
diffeomorphism. The domain $\mathcal{U}_{i}$ is called the coordinate domain
for the patch $\mathcal{P}^{0}_{i}$. Thus, the set
$\mathcal{A}^{0}=\\{(\mathcal{U}_{i},\eta_{i}^{0})\\}_{i=1}^{n_{p}}$ forms an
atlas for the manifold $\mathbb{S}^{2}$ [10] . We also know such an atlas
admits a smooth partition of unity subordinate to the open cover
$\\{\mathcal{P}_{i}^{0}\\}_{i=1}^{n_{p}}$ [37]. We choose a smooth partition
of unity $\\{\psi_{i}^{0}\\}_{i=1}^{n_{p}}$ where
$\psi_{i}^{0}:\mathbb{S}^{2}\longrightarrow\mathbb{R}$ such that
$\mathsf{supp}(\psi_{i}^{0})\subset\mathcal{P}_{i}^{0}$ for
$i=1,\ldots,n_{p}$, where $\mathsf{supp}(\cdot)$ denotes the support of the
function. Thus, for every $\bm{x}_{0}\in\mathbb{S}^{2}$,
$\sum_{i=1}^{n_{p}}\psi_{i}^{0}(\bm{x}_{0})=1$. The precise definition of
$\psi_{i}^{0}$ is given in Equation 3.9.

Define $\mathcal{P}_{i}:=\phi(\mathcal{P}_{i}^{0})$. Now
$\\{\mathcal{P}_{i}\\}_{i=1}^{n_{p}}$ form a set of overlapping patches that
cover the surface $\gamma$, _i.e.,_
$\bigcup_{i=1}^{n_{p}}\mathcal{P}_{i}=\gamma$. We then define coordinate maps
$\eta_{i}:\mathcal{U}_{i}\longrightarrow\mathcal{P}_{i}$, where
$\eta_{i}\equiv\phi\circ\eta_{i}^{0}$ and therefore, a corresponding atlas
$\mathcal{A}:=\\{(\mathcal{U}_{i},\eta_{i})\\}_{i=1}^{n_{p}}$ for the capsule
surface $\gamma$. A representation of the parameterization is shown in Figure
2. Since the patches are overlapping, a point $\bm{x}\in\gamma$ can belong to
multiple $\mathcal{P}_{i}$. For $\bm{x}\in\gamma$, let us define the set
$\mathcal{I}_{\bm{x}}=\\{1\leq i\leq n_{p}|\textit{
}\bm{x}\in\mathcal{P}_{i}\\}$. Thus, $\mathcal{I}_{\bm{x}}$ is the collection
of indices $i$ for which patch $\mathcal{P}_{i}$ contains $\bm{x}$. We also
define $\mathcal{P}_{ij}=\mathcal{P}_{i}\cap\mathcal{P}_{j}$ and
$\mathcal{U}_{ij}=\eta_{i}^{-1}(\mathcal{P}_{ij})$ for $i,j=1,\ldots,n_{p}$.
We also define transition maps
$\tau_{ij}:\mathcal{U}_{ij}\longrightarrow\mathcal{U}_{j}$ as
$\tau_{ij}\equiv(\eta_{j}^{0})^{-1}\circ\eta_{i}^{0}|_{\mathcal{U}_{ij}}$.
Furthermore, the diffeomorphism $\phi$ allows us to create a partition of
unity $\\{\psi_{i}\\}_{i=1}^{n_{p}}$ on $\gamma$ by defining it as

$\displaystyle\psi_{i}(\bm{x})=\psi_{i}^{0}(\phi^{-1}(\bm{x})),\quad\forall\bm{x}\in\gamma.$
(3.1)

For any function $f:\gamma\longrightarrow\mathbb{R}^{d}$, we can write

$\displaystyle
f(\bm{x})=\sum_{i=1}^{n_{p}}f(\bm{x})\psi_{i}(\bm{x}),\quad\forall\bm{x}\in\gamma.$
(3.2)

If $f$ is smooth, then $f\psi_{i}$ is smooth and compactly supported in
$\mathcal{P}_{i}$ for $i=1,\ldots,n_{p}$. We use this partition of unity
representation to compute derivatives and integrals on $\gamma$.

In this paper, we do not consider the adaptive case, and we assume that the
patch parameterization remains unchanged through the calculation. In our
numerical experiments, we used just six patches. Precisely, consider
$\mathcal{U}_{i}=(0,\pi)\times(0,\pi),\textit{ }i=1,\ldots,6$. The coordinate
maps $\\{\eta_{i}^{0}\\}_{i=1}^{6}$ are given below:

$\displaystyle\eta_{1}^{0}(u,v)$
$\displaystyle=(\sin{u}\cos{v},\sin{u}\sin{v},\cos{u}),$ (3.3)
$\displaystyle\eta_{2}^{0}(u,v)$
$\displaystyle=(-\sin{u}\cos{v},-\sin{u}\sin{v},\cos{u}),$ (3.4)
$\displaystyle\eta_{3}^{0}(u,v)$
$\displaystyle=(\sin{u}\sin{v},-\sin{u}\cos{v},\cos{u}),$ (3.5)
$\displaystyle\eta_{4}^{0}(u,v)$
$\displaystyle=(-\sin{u}\sin{v},\sin{u}\cos{v},\cos{u}),$ (3.6)
$\displaystyle\eta_{5}^{0}(u,v)$
$\displaystyle=(\sin{u}\cos{v},-\cos{u},\sin{u}\sin{v}),$ (3.7)
$\displaystyle\eta_{6}^{0}(u,v)$
$\displaystyle=(\sin{u}\cos{v},\cos{u},-\sin{u}\sin{v}).$ (3.8)

These six charts form six hemispherical patches
$\\{\mathcal{P}_{i}^{0}\\}_{i=0}^{6}$ that cover the unit sphere as shown in
Figure 3. The computation of the transition maps $\tau_{ij}$ is relatively
simple and is given in the Section A.1. We use the bump function,
$b(r)=\frac{e^{2e^{-1/|r|}}}{|r|-1}$ for $|r|<1$ and 0 otherwise, to construct
the partition of unity functions on each patch. For the particular
parameterization we use, for each patch of the unit sphere, we define

$\displaystyle\psi_{i}^{0}(\bm{x}_{0})=\frac{b\left(\frac{d(\bm{x}_{0},\eta_{i}^{0}(\pi/2,\pi/2))}{r_{0}}\right)}{\sum_{j=1}^{n_{p}}b\left(\frac{d(\bm{x}_{0},\eta_{j}^{0}(\pi/2,\pi/2))}{r_{0}}\right)},\
\forall\bm{x}_{0}\in\mathbb{S}^{2}$ (3.9)

where $d(\bm{x}_{0},\bm{y}_{0})$ denotes the great circle distance between
$\bm{x}_{0}$ and $\bm{y}_{0}$ on the unit sphere, and $r_{0}>0$ determines the
support of the partition of unity inside the patch. The argument of $b(\cdot)$
is the normalized great circle distance from a point $\bm{x}$ to the center of
the patch $\eta_{i}^{0}(\pi/2,\pi/2)$, where the normalizing factor $r_{0}$ is
chosen to be $r_{0}=5\pi/12$ (see Section A.3). Now, for a given surface
$\gamma$ and diffeomorphism $\phi$, we readily have an atlas and the
transition maps for the parameterization of $\gamma$. The corresponding
partition of unity on $\gamma$ is available using Equation 3.1 and Equation
3.9.

Figure 2: Here we summarize the notation for the atlas construction. The unit
sphere $\mathbb{S}^{2}$ is on the left and the capsule $\gamma$ is on the
right with the diffeomorphism $\phi:\mathbb{S}^{2}\longrightarrow\gamma$. We
show two overlapping patches colored in red and blue. The first patch
$\mathcal{P}_{1}^{0}$ is the red colored arc from point $a_{0}$ to $b_{0}$ on
$\mathbb{S}^{2}$. Its corresponding patch $\mathcal{P}_{1}$ on the capsule
surface $\gamma$ is shown in red as the arc from points $a_{1}$ to $b_{1}$.
The second patch $\mathcal{P}_{2}^{0}$ is the blue colored arc from point
$c_{0}$ to $d_{0}$ on $\mathbb{S}^{2}$. Its corresponding patch
$\mathcal{P}_{2}$ on the capsule surface $\gamma$ is shown in blue as the arc
from points $c_{1}$ to $d_{1}$. Their corresponding coordinate domains
$\mathcal{U}_{1}$ and $\mathcal{U}_{2}$ are also shown in the red and blue
color. The coordinate charts
$\eta_{i}^{0}:\mathcal{U}_{1}\longrightarrow\mathcal{P}^{0}_{i},\ i=1,2$, are
shown as dashed lines in the respective colors. The diffeomorphism $\phi$ also
gives the coordinate charts
$\eta_{i}:\mathcal{U}_{1}\longrightarrow\mathcal{P}_{i},\ i=1,2$, for the
patches on the capsule surface $\gamma$ shown as colored dashed lines from
$\mathcal{U}_{i}$ to $\mathcal{P}_{i}$. The corresponding partition of unity
functions $\psi^{0}_{i},i=1,2,$ with
$\mathsf{supp}(\psi^{0}_{i})\subset\mathcal{P}^{0}_{i}$ are drawn over
coordinate domains $\mathcal{U}_{i}$ for visual clarity since
$\mathcal{P}_{i}^{0}$ is diffeomorphic to $\mathcal{U}_{i}$.

Figure 3: The six hemispherical patches forming an open cover of the unit
sphere $\mathcal{S}^{2}$. Each one is represented by the black grid. a)–f)
$\mathcal{P}_{i}^{0}$ for $i=1,2,\ldots,6$.

###  3.2 Surface discretization

We now describe in detail the discretization of the surface $\gamma$ using the
parameterization described above in Section 3.1. As mentioned above,
$\mathcal{U}_{i}=(0,\pi)\times(0,\pi)$ for our parameterization. We use the
following $m_{\mathrm{th}}$-order uniform grid in $\mathcal{U}_{i}$ as
follows:

$\displaystyle U^{m,i}_{j,k}=\left(\frac{j\pi}{m},\frac{k\pi}{m}\right),\
\forall j,k\in\\{1,\ldots,m-1\\}.$ (3.10)

We use $h_{m}$ to denote the $m_{\mathrm{th}}$-order grid spacing in each
$\mathcal{U}_{i}$ space, _i.e.,_ $h_{m}=\frac{\pi}{m}$. The discretization
points for each patch $\mathcal{P}_{i}^{0}$ on the unit sphere are given by

$\displaystyle X^{0,m,i}_{j,k}=\eta_{i}^{0}(U^{m,i}_{j,k}),\ \forall
j,k\in\\{1,\ldots,m-1\\}.$ (3.11)

The discretization points for each patch $\mathcal{P}_{i}$ on $\gamma$ are
given by

$\displaystyle X^{m,i}_{j,k}=\eta_{i}(U^{m,i}_{j,k}),\ \forall
j,k\in\\{1,\ldots,m-1\\}.$ (3.12)

Thus, each patch contains $(m-1)^{2}$ discretization points for an
$m_{\mathrm{th}}$-order grid leading to a total of $N=6(m-1)^{2}$ points for
the surface. The dynamics of the capsule are represented as the time
trajectories $X^{m,i}_{j,k}(t)$. The partition of unity values at the
discretization points $\psi^{m,i}_{j,k}=\psi_{i}(X^{m,i}_{j,k})$ can be
precomputed and stored to be used later for computing integrals and
derivatives (see Equation 3.1). The sample grids for a unit sphere and an
ellipsoid are given in Figure 4 and Figure 5. To get the diffeomorphism $\phi$
for an ellipsoid $\gamma$ given by the level surface
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, we consider
the spherical angles parameterization
$\beta(u,v):[0,\pi]\times[0,2\pi)\longrightarrow\gamma$ given by

$\displaystyle\beta(u,v)=(a\sin{u}\cos{v},b\sin{u}\sin{v},c\cos{u}),u\in[0,\pi],v\in[0,2\pi).$
(3.13)

We also consider the spherical angles parameterization of the unit sphere
$\beta_{0}(u,v):[0,\pi]\times[0,2\pi)\longrightarrow\mathbb{S}^{2}$ given by

$\displaystyle\beta_{0}(u,v)=(\sin{u}\cos{v},\sin{u}\sin{v},\cos{u}),u\in[0,\pi],v\in[0,2\pi).$
(3.14)

The diffeomorphism $\phi:\mathbb{S}^{2}\longrightarrow\gamma$ is then given as

$\phi(\bm{x}_{0})=\begin{cases}\beta(\beta_{0}^{-1}(\bm{x}_{0}))&\text{ if
}\bm{x}_{0}\in\mathbb{S}^{2}\backslash\\{(0,0,-1),(0,0,1)\\},\\\
(0,0,-c)&\text{ if }\bm{x}_{0}=(0,0,-1),\\\ (0,0,c)&\text{ if
}\bm{x}_{0}=(0,0,1).\\\ \end{cases}$

_Remark 1:_ Given an initial capsule surface $\gamma$, we only use the
diffeomorphism $\phi:\mathbb{S}^{2}\longrightarrow\gamma$ to compute the atlas
$\mathcal{A}$, the partition of unity $\\{\psi_{i}\\}_{i=1}^{n_{p}}$ for
$\gamma$ and the transition maps $\tau_{ij}$ for the coordinate domains. Once
computed, we initialize the surface discretization points (see Equation 3.12),
and track their time trajectories to simulate the dynamics of the capsule. We
do not need to use $\phi$ afterwards for capsule dynamics simulation. We do
not change the partition of unity values after initialization. Hence, we have
$\psi^{m,i}_{j,k}=\psi_{i}(X^{m,i}_{j,k}(t_{0}))=\psi_{i}(X^{m,i}_{j,k}(t))$
where $t_{0}$ is the initial time. These partition of unity values are
computed at initialization and stored for subsequent use. Also, the transition
maps between the coordinate domains remain unchanged and are precomputed for
usage later.

_Notation:_ For brevity, we use $\bm{U}^{m}$ to denote the $N\times 2$ matrix
containing all the $m_{\mathrm{th}}$-order grid points, _i.e.,_
$\bm{U}^{m}=[\\{U^{m,i}_{j,k}\\}_{1\leq j,k\leq m,1\leq i\leq n_{p}}]^{T}$.
Similarly, we define $\bm{X}^{m}=[\\{X^{m,i}_{j,k}\\}_{1\leq j,k\leq m,1\leq
i\leq n_{p}}]^{T}$ to be the $N\times 3$ matrix of corresponding
discretization points on $\gamma$,
$\bm{\Psi}^{m}=[\\{\psi^{m,i}_{j,k}\\}_{1\leq j,k\leq m,1\leq i\leq
n_{p}}]^{T}$ to be the $N\times 1$ column vector containing the values of
partition of unity functions and $\bm{F}^{m}$ to be the $N\times 3$ matrix
containing the interfacial force $\bm{f}$ values at discretization points of
the $m_{\mathrm{th}}$-order grid. Further, we use $\bm{U}^{m,i}$ to denote the
$(m-1)^{2}\times 2$ matrix of grid points belonging to $\mathcal{U}_{i}$.
Similarly, we use $\bm{X}^{m,i},\bm{\Psi}^{m,i},\bm{F}^{m,i}$ to denote the
values in the patch $\mathcal{P}_{i}$. For the unit sphere, we denote the
corresponding discretization points as $\bm{X}^{0,m}$. The partition of unity
column vector on $\bm{X}^{0,m}$ is the same as $\bm{\Psi}^{m}$.

_Remark 2:_ While we do not need $\phi$ for simulating capsule dynamics after
the initialization of atlas, transition maps and partition of unity, we use
$||\nabla_{\mathbb{S}^{2}}\phi||_{\infty}$ as a quality metric to gauge the
smoothness of the capsule surface $\gamma$ during the simulations in this
paper (see Section 3.6). To compute
$||\nabla_{\mathbb{S}^{2}}\phi||_{\infty}$, we will need the discretization
points $\bm{X}^{0,m}$. Hence, we store these values as well. Note that
$\phi(\bm{X}^{0,m})=\bm{X}^{m}$.

(a) $m=8$

(b) $m=16$

(c) $m=32$

Figure 4: Representation of discretization of a unit sphere using
$m_{\mathrm{th}}$-order grids for (a) $m=8$, (b) $m=16$, (c) $m=32$.

(a) $m=8$

(b) $m=16$

(c) $m=32$

Figure 5: Representation of discretization of the ellipsoid
$x^{2}/a^{2}+y^{2}/b^{2}+z^{2}/c^{2}=1$ with $a=0.5,b=1,c=1$, using
$m_{\mathrm{th}}$-order grids for (a) $m=8$, (b) $m=16$, (c) $m=32$.

###  3.3 Smooth surface integrals

Let $f:\gamma\longrightarrow\mathbb{R}^{d}$ be a smooth function. We can write
its surface integral as a sum of the integral over all the patches using the
partition of unity $\\{\psi_{i}\\}_{i=1}^{n_{p}}$. Then, we use the product
periodic trapezoidal rule using the function values at the grid points
mentioned in Equation 3.12, which gives superalgebraic convergence [7]. We
write

$\displaystyle\int_{\gamma}fd\gamma=\sum_{i=1}^{n_{p}}\int_{\mathcal{P}_{i}}f\psi_{i}d\gamma=\sum_{i=1}^{n_{p}}\int_{\mathcal{U}_{i}}f\psi_{i}W_{i}d\mathcal{U}_{i}\approx\sum_{i=1}^{n_{p}}\left(\sum_{j,k\in\\{1,\ldots,m-1\\}}f(X^{m,i}_{j,k})\psi^{m,i}_{j,k}W_{i}(X_{j,k}^{m,i})h_{m}^{2}\right),$
(3.15)

where $W_{i}$ is the surface area element for patch $\mathcal{P}_{i}$
(numerical computation of $W_{i}$ is discussed in Section 3.5).

_Work complexity:_ The work complexity to compute this integral numerically
for $m_{\mathrm{th}}$-order grid is $O(n_{p}m^{2}).$ In our specific
parameterization, $n_{p}=6$. Hence, the work required is $O(6m^{2})$.

###  3.4 Singular integration

We use the regularized Stokes kernel described in [32] to evaluate Stokes
potentials to high-order accuracy using the discretization described above in
Section 3.2.

Consider a smooth vector-valued function
$\bm{f}:\gamma\longrightarrow\mathbb{R}^{3}$. The single-layer Stokes
potential of $\bm{f}$ at $\bm{x}\in\gamma$ as described in Section 2.2 can be
written as

$\displaystyle\mathcal{S}_{\gamma}[\bm{f}](\bm{x})=\frac{1}{8\pi\mu}\int_{\gamma}\left(\frac{\bm{f}(\bm{y})}{r}+(\bm{f}(\bm{y})\cdot(\bm{x}-\bm{y}))(\bm{x}-\bm{y})\frac{1}{r^{3}}\right)d\gamma,$
(3.16)

where $r=||\bm{x}-\bm{y}||$. The regularized version of Stokes single layer
potential [32] is given as

$\displaystyle\mathcal{S}_{\gamma}^{\delta}[\bm{f}](\bm{x})=\frac{1}{8\pi\mu}\int_{\gamma}\left(\frac{\bm{f}(\bm{y})s_{1}(r/\delta)}{r}+(\bm{f}(\bm{y})\cdot(\bm{x}-\bm{y}))(\bm{x}-\bm{y})\frac{s_{2}(r/\delta)}{r^{3}}\right)d\gamma,$
(3.17)

where $\delta>0$ is the regularization parameter and $s_{1},s_{2}$ are
smoothing factors such that $\frac{s_{1}(r/\delta)}{r}$ and
$\frac{s_{2}(r/\delta)}{r^{3}}$ are smooth functions as $r\longrightarrow 0$.
These smoothing factors ensure that $S_{\gamma}^{\delta}[\bm{f}]$ is an
integral involving a smooth kernel and can be evaluated using the periodic
trapezoidal rule as in Section 3.3. Following [32, 24], we choose the
smoothing factors $s_{1}$ and $s_{2}$ to be

$\displaystyle s_{1}(r)$
$\displaystyle=\operatorname{Erf}(r)-\frac{2}{3}r(2r^{2}-5)\frac{e^{-r^{2}}}{\sqrt{\pi}},$
(3.18) $\displaystyle s_{2}(r)$
$\displaystyle=\operatorname{Erf}(r)-\frac{2}{3}r(4r^{4}-14r^{2}+3)\frac{e^{-r^{2}}}{\sqrt{\pi}},$
(3.19)

where $\operatorname{Erf}(\cdot)$ is the error function. Using Taylor
expansion, we can show that as $r\longrightarrow 0$,
$\frac{s_{1}(r/\delta)}{r}\longrightarrow\frac{16}{3\delta\sqrt{\pi}}$ and
$\frac{s_{2}(r/\delta)}{r^{3}}\longrightarrow\frac{32}{3\delta^{3}\sqrt{\pi}}$.
This choice of smoothing factors ensures that the regularized Stokes potential
is $O(\delta^{5})$ accurate [32]. This regularized integral is discretized by
the trapezoidal rule for smooth functions. However, the question is what
should be the relation between $m$ and $\delta$?

_Regularization parameter $\delta$:_ The choice of regularization parameter
$\delta$ is crucial. As discussed in [4, 32], the chosen $\delta$ should be
large enough that regularization error dominates the discretization error in
computing the integral in Equation 3.17. In our simulations, when computing
the Stokes potential for a target point $X^{m,i}_{j,k}$ in patch
$\mathcal{P}_{i}$, we choose a $\mathcal{P}_{i}$-dependent regularization
parameter $\delta_{i}$ at every time step as the capsule evolves. It is
defined by

$\displaystyle\delta^{m,i}_{j,k}=C\delta^{*}\mbox{~{}where~{}}\delta^{*}=\left(\max_{l,l^{\prime}\in\\{1,\ldots,m-1\\}}\\{d_{e}(X^{m,i}_{l,l^{\prime}},\mathcal{X}(X^{m,i}_{l,l^{\prime}})\\}\right).$
(3.20)

Here, we choose the constant $C=1$ (see Section A.7 for more details on the
choice $C$),
$\mathcal{X}(X^{m,i}_{l,l^{\prime}})=\\{X^{m,i}_{l-a,l^{\prime}-b}:a,b\in\\{-1,0,1\\}\\}$
refers to the set of points which are adjacent neighbors of
$X^{m,i}_{l,l^{\prime}}$ and $d_{e}(a,B)$ denotes the maximum Euclidean
distance between $a$ and the points in set $B$. For brevity, we define
$\bm{\delta}^{m}=[\\{\delta^{m,i}_{j,k}\\}_{1\leq j,k\leq m-1,1\leq i\leq
n_{p}}]^{T}$. We experimentally observed that our choice for $\delta$ improves
the accuracy the singular quadrature as in our simulations as capsule changes
shapes and patches evolve over time, as opposed to using a fixed global
regularization parameter $\delta=C^{\prime}h_{m}$ where $C^{\prime}$ is a
positive constant (like in [32]). We tabulate the relative errors for
different values of $C^{\prime}$ in using fixed global regularization
parameter $\delta=C^{\prime}h_{m}$ on a sample ellipsoidal shape in Table
11(a) to illustrate this.

_Upsampling:_ In our experiments, the observed prefactor in our weakly
singular quadrature scheme is quite significant. For example at small $m$, say
$m=16$, we can resolve the surface well but we need more quadrature points to
resolve the weakly singular integrals. For this reason we use upsampling,
especially in the small-$m$ regime. Using numerical experiments, we determined
that a four times upsampling works well. Thus, we use an upsampled grid with
$N_{\mathrm{up}}=6(4m-1)^{2}$ grid points for evaluating the single layer
potential. We use cubic spline interpolation [9] on each patch to upsample the
surface discretization points $\bm{X}^{m}$ and surface interfacial force
$\bm{F}^{m}$. We define this upsampling operator as $\bm{I}_{u}^{m}$. Thus,
the upsampled surface points and surface force vector can be written as

$\displaystyle\bm{X}^{m}_{u}=\bm{I}_{u}^{m}\bm{X}^{m},\bm{F}^{m}_{u}=\bm{I}_{u}^{m}\bm{F}^{m}.$
(3.21)

The partition of unity values on the upsampled grid, denoted by
$\bm{\Psi}^{m}_{u}$, can also be precomputed and stored for use. We use the
upsampled values to compute the Stokes single layer potential ( Equation 3.17)
and then downsample the Stokes potential to get the values on
$m_{\mathrm{th}}$-order grid. Note that both $\bm{I}_{u}^{m}$ and
$\bm{I}_{d}^{m}$ are cubic spline based weighted interpolation operators.
Their matrices can be precomputed and stored for use throughout the
simulation. We also compute the $\bm{\delta}^{m}_{u}$ for the upsampled grid
using Equation 3.20.

_Work complexity:_ The work complexity for the evaluation of Stokes single
layer potential on upsampled grid for a target point is $O(96m^{2})$. For
$N_{up}$ target points in upsampled grid, the total work complexity for
computing single layer potential is $O(9216m^{4})$. The complexity for
upsampling and downsampling is $O(m^{2})$ since $\bm{I}_{u}^{m}$ and
$\bm{I}_{d}^{m}$ are sparse linear operators. Hence, the total work complexity
is $O(m^{4})$.

_GPU acceleration:_ The evaluation of single layer potential on upsampled grid
is the most expensive part of our numerical scheme. In order to accelerate
this computation, we use a CUDA implementation of computation of Stokes single
layer potential computation based on the all-pairs $O(N^{2})$ calculation
using CUDA on GPU (for details refer to [25]).

###  3.5 Surface derivatives

In this section we discuss in detail our numerical scheme for calculating
surface derivatives. The derivatives are needed for computing the interfacial
force $\bm{f}$ and the surface area element $W$ (for Equation 3.15) at each
discretization point. The first step in the computation of the interfacial
force $\bm{f}$ is the computation of the shear stress tensor $\Lambda$. The
computation of $\Lambda$ requires the surface tangents and normals in the
current configuration (see Equation 2.8). After computing the stress tensor
$\Lambda$, the second step is to compute the interfacial force by calculating
the surface divergence of $\Lambda$ as $\bm{f}=\nabla_{\gamma}\cdot\Lambda$.
For the computation of surface area element $W$, we need the coefficients
$E,F,G$ of the first fundamental form of the surface. The formulas for the
surface divergence, denoted by $\nabla_{\gamma}\cdot$, and for the surface
area element $W$ along with $E,F,G$ are summarized in the Section A.2.

Let $\eta(u,v):\mathcal{U}\longrightarrow\mathcal{P}\subset\gamma$ be a
surface patch of $\gamma$ where $\eta\in\\{\eta_{1},\ldots,\eta_{n_{p}}\\}$.
Let $\eta(u_{0},v_{0})=\bm{x}\in\mathcal{P}$. Below, we summarize the steps
required to evaluate the required shape derivatives
$\bm{f}(\eta(u_{0},v_{0}))$ and $W(\eta(u_{0},v_{0}))$.

  1. 1.

We need to compute the tangents, denoted by $\bm{x}_{u}$ and $\bm{x}_{v}$, and
the unit normal $\bm{n}(\bm{x})$ for the current and the reference
configuration. . These quantities are given by
$\bm{x}_{u}=\frac{\partial\eta}{\partial
u}|_{(u_{0},v_{0})},\bm{x}_{v}=\frac{\partial\eta}{\partial
v}|_{(u_{0},v_{0})}$ and
$\bm{n}(\bm{x})=\frac{\bm{x}_{u}\times\bm{x}_{v}}{||\bm{x}_{u}\times\bm{x}_{v}||}$.
Given these derivatives we can then compute $\Lambda(\eta(u_{0},v_{0}))$.

  2. 2.

The next step is to compute the surface divergence of $\Lambda(\eta(u,v))$ at
$(u_{0},v_{0})$. This involves computation of $E,F,G$ and partial $u,v$
derivatives of the components of $\Lambda(\eta(u,v))$ (see Section A.2). The
computation of $E,F,G$ at point $\bm{x}$ is straightforward using $\bm{x}_{u}$
and $\bm{x}_{v}$. The surface area element is then computed using
$W(\bm{x}(u_{0},v_{0}))=\sqrt{EG-F^{2}}$.

We note here from the above summary that the computation of both $\bm{f}$ and
$W$ boils down to multiple computations of the $u$ and $v$ partial derivatives
of functions like the coordinate chart $\eta(u,v)$ and the components of
$\Lambda(\eta(u,v))$. Now given an arbitrary function
$g:\gamma\longrightarrow\mathbb{R}$, we consider the parametric function
$\tilde{g}^{i}\equiv
g(\eta_{i}(u,v)):\mathcal{U}_{i}\longrightarrow\mathbb{R}$ for
$i=1,\ldots,n_{p}$. We denote its partial $u,v$ derivatives as
$\tilde{g}^{i}_{u}$ and $\tilde{g}^{i}_{v}$. We describe below the numerical
scheme to compute these partial derivatives $\tilde{g}^{i}_{u}$ and
$\tilde{g}^{i}_{v}$. This scheme is then sufficient to compute $\bm{f}$ and
$W$ at the discretization grid $\bm{U}^{m}$ using the formulas mentioned in
Section A.2. Before we describe the scheme in detail, we summarize the three
main steps of the scheme below:

Figure 6: A 2D representation of a patch extension to apply the central finite
difference stencil. The unit sphere $\mathbb{S}^{2}$ is on the left and the
capsule $\gamma$ is on the right with the diffeomorphism
$\phi:\mathbb{S}^{2}\longrightarrow\gamma$. The discrete points on the patch
$\mathcal{P}_{2}^{0}$ are represented by the blue colored dotted arc from
point $c_{0}$ to $d_{0}$ on $\mathbb{S}^{2}$. Its corresponding patch
$\mathcal{P}_{2}$ on the capsule surface $\gamma$ is shown in dotted blue line
as the arc from points $c_{1}$ to $d_{1}$. Its corresponding coordinate domain
$\mathcal{U}_{2}$ is also shown in the blue color. The extended grid points
are shown in ’$+$’ symbol alongside $\mathcal{U}_{2}$. The extended
discretization points are on $\mathbb{S}^{2}$ and $\gamma$ are shown with the
’$\times$’ symbol. The coordinate map $\eta_{2}^{0}$ we use in our
parameterization has a natural extension to the extended domain with extended
grid points with the same expression as in Equation 3.8.

  1. 1.

_Patch extension:_ We intend to use central finite difference stencil to
compute $\tilde{g}^{i}_{u}$ and $\tilde{g}^{i}_{v}$ at the discretization
points for $i=1,\ldots,n_{p}$. To apply central difference at the grid points
near the boundary of each coordinate domain $\mathcal{U}_{i}$, we need
function values on ghost grid points outside the discretization grid in
$\mathcal{U}_{i}$. We call this process of obtaining values on the extended
grid points as _patch extension_. We obtain the values of $\tilde{g}^{i}$ on
these ghost grid points using the patch extension process. Let
$\bm{x}_{1}=\eta_{i}(u_{1}^{i},v_{1}^{i})$ for $i=1,\ldots,n_{p}$, be an
extended discretization point (see Figure 6). The value of the function
$\tilde{g}^{i}$ at the extended grid point $(u_{1}^{i},v_{1}^{i})$, is given
in the continuous form as the weighted average of the values across the
patches as follows:

$\displaystyle\tilde{g}^{i}(u_{1}^{i},v_{1}^{i})=\sum_{1\leq j\leq
n_{p}}\tilde{g}^{j}(u_{1}^{j},v_{1}^{j})\psi_{j}(\bm{x}_{1}).$ (3.22)

Here, the partition of unity values are used as the weights for each patch. In
the continuous case,
$\tilde{g}^{i}(u_{1}^{i},v_{1}^{i})=\tilde{g}^{j}(u_{1}^{j},v_{1}^{j})$ for
$i,j=1,\ldots,n_{p}$. Hence, it is easy to see that the Equation 3.22 is valid
in the continuous setting by putting
$\tilde{g}^{i}(u_{1}^{i},v_{1}^{i})=\tilde{g}^{j}(u_{1}^{j},v_{1}^{j})$ for
$i,j=1,\ldots,n_{p}$ $i,j=1,\ldots,n_{p}$ and noting that the partition of
unity values add up to unity. In the discrete case, the above equation
provides a unique definition of discretized function values at every extended
grid point since $\tilde{g}^{i}(u_{1}^{i},v_{1}^{i})$ and
$\tilde{g}^{j}(u_{1}^{j},v_{1}^{j})$ are not necessarily the same for
$i,j\in\\{1,\ldots,n_{p}\\}$ in the discrete case due to numerical errors.

  2. 2.

_Finite difference:_ After extension of function values to ghost grid points,
we use central finite difference operator on extended grids to obtain
$\tilde{g}^{i}_{u}$ and $\tilde{g}^{i}_{v}$ on the grid points $\bm{U}^{m,i}$
for $i=1,\ldots,n_{p}$.

  3. 3.

_Blending:_ We note that a point $\bm{x}\in\gamma$ can belong to several
different patches. We suppose $\bm{x}=\eta_{i}(u_{0}^{i},v_{0}^{i})$ for
$i=1,\ldots,n_{p}$. The interfacial force $\bm{f}(\bm{x})$ is intrinsic to the
surface $\gamma$ and is independent of the local surface parameterization.
Using patch parameterization dependent numerical derivatives
$\tilde{g}^{i}_{u}$ and $\tilde{g}^{i}_{v}$ leads to different discrete
derivative values at $(u_{0}^{i},v_{0}^{i})$ for different $i=1,\ldots,n_{p}$.
To get unique derivative values consistent across all the patches, we _blend_
the derivatives across all the patches. The blending process is taking the
weighted average of derivative values across all the patches except that we
also have to transform the derivatives from, say the domain $\mathcal{U}_{j}$
to the domain $\mathcal{U}_{i}$. This transformation is given by the Jacobian
$\bm{J}_{\tau_{ij}}$ of the transition map
$\tau_{ij}(u,v)=(\tau_{ij}^{(1)}(u,v),\tau_{ij}^{(2)}(u,v))\in\mathcal{U}_{j}$
as follows:

$\displaystyle\begin{bmatrix}\tilde{g}_{u}^{i}\\\
\tilde{g}_{v}^{i}\end{bmatrix}=\bm{J}_{\tau_{ij}}\begin{bmatrix}\tilde{g}_{u}^{j}\\\
\tilde{g}_{v}^{j}\end{bmatrix},$ (3.23)

where the Jacobian matrix is

$\displaystyle\bm{J}_{\tau_{ij}}=\begin{bmatrix}\partial\tau_{ij}^{(1)}/\partial
u&\partial\tau^{(2)}_{ij}/\partial u\\\ \partial\tau_{ij}^{(1)}/\partial
v&\partial\tau_{ij}^{(2)}/\partial v\end{bmatrix}.$ (3.24)

Here, $\tau_{ij}^{(1)},\tau_{ij}^{(2)}$ are defined as the first and second
component respectively of the transition map, _i.e.,_
$\tau_{ij}(u,v)=(\tau_{ij}^{(1)}(u,v),\tau_{ij}^{(2)}(u,v))\in\mathcal{U}_{j}$.
Both the RHS and the LHS in Equation 3.23 are evaluated at
$(u^{i}_{0},v^{i}_{0})$. Using this Jacobian transformation, we define the
blending process as the weighted average using partition of unity values as
weights for the patches:

$\displaystyle\begin{bmatrix}\tilde{g}_{u}^{i}\\\
\tilde{g}_{v}^{i}\end{bmatrix}=\sum_{j=1}^{n_{p}}\psi_{j}(\bm{x})\bm{J}_{\tau_{ij}}\begin{bmatrix}\tilde{g}_{u}^{j}\\\
\tilde{g}_{v}^{j}\end{bmatrix}$ (3.25)

for all $i=1,\ldots,n_{p}$, where $\bm{x}=\eta_{i}(u^{i}_{0},v^{i}_{0})$. The
partition of unity values $\psi_{j}$ form the respective weights for the
coordinate domain $\mathcal{U}_{i}$ like in the Equation 3.22. Notice that the
Equation 3.25 is trivially valid in the continuous case since Equation 3.23
holds and $\psi_{j}$ sum to unity. The above equation can be expanded as
follows:

$\displaystyle\tilde{g}^{i}_{u}(u^{i}_{0},v^{i}_{0})=\tilde{g}^{i}_{u}(u^{i}_{0},v^{i}_{0})\psi_{i}(\bm{x})+\sum_{j\neq
i,1\leq j\leq
n_{p}}\left(\tilde{g}_{u}^{j}\left.\frac{\partial\tau_{ij}^{(1)}}{\partial
u}\right|_{(u^{i}_{0},v^{i}_{0})}+\tilde{g}_{v}^{j}\left.\frac{\partial\tau_{ij}^{(2)}}{\partial
u}\right\rvert_{(u^{i}_{0},v^{i}_{0})}\right)\psi_{j}(\bm{x})$ (3.26)
$\displaystyle\tilde{g}^{i}_{v}(u^{i}_{0},v^{i}_{0})=\tilde{g}^{i}_{v}(u^{i}_{0},v^{i}_{0})\psi_{i}(\bm{x})+\sum_{j\neq
i,1\leq j\leq
n_{p}}\left(\tilde{g}_{u}^{j}\left.\frac{\partial\tau_{ij}^{(1)}}{\partial
v}\right\rvert_{(u^{i}_{0},v^{i}_{0})}+\tilde{g}_{v}^{j}\left.\frac{\partial\tau_{ij}^{(2)}}{\partial
v}\right\rvert_{(u^{i}_{0},v^{i}_{0})}\right)\,\psi_{j}(\bm{x}).$ (3.27)

The terms inside the parenthesis in Equation 3.26 and Equation 3.27 come from
the Jacobian transformation in Equation 3.23. We list all the transition maps
for the specific parameterization we use in Section A.1. Equation 3.26 and
Equation 3.27 describe the blending equation we use to get the weighted
average of derivatives across all patches.

We now describe in detail all the three steps mentioned above and write the
discretized versions of these steps.

_Notation:_ We use $\tilde{\bm{g}}^{m}$ to denote the column vector of the
values of $\tilde{g}^{i}$ at the grid points $\bm{U}^{m,i}$ for all
$i=1,\ldots,n_{p}$. We use $\tilde{\bm{g}}^{m,i}$ to denote the column vector
of the values at the grid points $\bm{U}^{m,i}$.

_Patch extension:_ As discussed above, for the discretization points near the
boundary of $\mathcal{U}_{i},i=1,\ldots,n_{p}$, we need function values on
_ghost points_ lying outside $\mathcal{U}_{i}$ in order to use the central FDM
stencil. In our simulation, we use a fifth order 7-point symmetric 1D-stencil
[16] for each of the partial derivatives. Thus, in order to calculate the
derivative values on the grid $\bm{U}^{m}$, we need the values of $\tilde{g}$
on the extended grid (including ghost nodes) given by

$\displaystyle U^{m,i,ext}_{j,k}=\left(\frac{j\pi}{m},\frac{k\pi}{m}\right),\
\forall j,k\in\\{-2,-1,\ldots,m+1,m+2\\},$ (3.28)

where the set of points defined by
$G^{m,i}:=\\{U^{m,i,ext}_{j,k}\\}_{j,k\in\\{-2,-1,\ldots,m+1,m+2\\}}\backslash\\{U^{m,i,ext}_{j,k}\\}_{j,k\in\\{1,\ldots,m-1\\}}$
denotes the ghost discretization points for $\mathcal{U}_{i}$. We define the
extended coordinate chart domain as
$\mathcal{U}^{m,ext}_{i}=(\frac{-3\pi}{m},\frac{(m+3)\pi}{m})\times(\frac{-3\pi}{m+1},\frac{(m+4)\pi}{m+1})$
containing the extended grid points as defined in Equation 3.28. The
parameterization of unit sphere $\eta_{i}^{0}$ has a natural extension to
$\mathcal{U}_{i}^{m,ext}$ with the same function expressions as given in
Equation 3.8. Using the diffeomorphism $\phi$, we also have mappings
$\eta_{i}^{m,ext}:\mathcal{U}_{i}^{m,ext}\longrightarrow\gamma$. We define
$\mathcal{P}_{i}^{m,ext}:=\eta_{i}^{m,ext}(\mathcal{U}_{i}^{m,ext})$. We
define the matrix of extended grid points for all the coordinate domains as
$\bm{U}^{m,ext}=[\\{U^{m,i,ext}_{j,k}\\}_{-2\leq j,k\leq m+2,1\leq i\leq
n_{p}}]^{T}$, while $\bm{U}^{m,i,ext}$ refers to the matrix of extended grid
points for the coordinate domain $\mathcal{U}_{i}$. The corresponding
discretization points on the surface $\gamma$ are theoretically given by

$\displaystyle X^{m,i,ext}_{j,k}=\eta_{i}\left(U^{m,i,ext}_{j,k}\right),\
\forall j,k\in\\{-2,\ldots,m+2\\}.$ (3.29)

Let the matrix of all surface discretization nodes on extended patch be
$\bm{X}^{m,ext}=[\\{X^{m,i}_{j,k}\\}_{-2\leq j,k\leq m+2,1\leq i\leq
n_{p}}]^{T}$. Our goal of patch extension is to find the values of the
function $\tilde{g}^{i}$ on $\bm{U}^{m,i,ext}$, denoted by
$\tilde{\bm{g}}^{m,i,ext}$, using the known values on $\bm{U}^{m,i}$. The
discretized form of patch extension is given as

$\displaystyle\tilde{\bm{g}}^{m,i,ext}$ $\displaystyle=\sum_{1\leq j\leq
n_{p}}\left(\bm{I}_{ij}^{m,ext}\tilde{\bm{g}}^{m,j}\right)\left.\psi_{j}\right\rvert_{\bm{X}^{m,i,ext}}.$
(3.30)

This equation is the discretized version of the Equation 3.22. Here, we
additionally use the cubic splines interpolation, denoted by the operator
$\bm{I}_{ij}^{m,ext}$, to get the values of $\tilde{g}^{j}$ on
$\bm{U}^{m,i,ext}$ using the known values on the discrete grid $\bm{U}^{m,j}$.

_Finite Differences:_ We use $\bm{D}^{m,i}_{u}$ and $\bm{D}^{m,i}_{v}$ to
denote the central finite difference operators for computing partial
derivative with respect to $u$ and $v$ respectively on the
$m_{\mathrm{th}}$-order grid in coordinate chart domain $\mathcal{U}_{i}$.
Thus, we write the partial derivatives on $\bm{U}^{m,i}$ as

$\displaystyle\tilde{\bm{g}}_{u}^{m,i}=\bm{D}^{m,i}_{u}\tilde{\bm{g}}^{m,i,ext},\tilde{\bm{g}}_{v}^{m,i}=\bm{D}^{m,i}_{v}\tilde{\bm{g}}^{m,i,ext},$
(3.31)

where $\tilde{\bm{g}}^{m,i,ext}$ is calculated using Equation 3.30.

_Blending:_ As discussed above, applying Equation 3.31 can lead to different
numerical values of derivatives at the same discretization point for different
domains $\mathcal{U}_{i}$. Therefore, we use blending to get unique values
across all the patches. The discretized versions of the blending equations
Equation 3.26 and Equation 3.27 are given as follows:

$\displaystyle\tilde{\bm{g}}_{u}^{m,i}=\tilde{\bm{g}}_{u}^{m,i}\left.\psi_{i}\right\rvert_{\bm{X}^{m,i}}+\sum_{j\neq
i,1\leq j\leq
n_{p}}\left(\left(\bm{I}_{ij}^{m}\tilde{\bm{g}}_{u}^{m,j}\right)\left.\frac{\partial\tau_{ij}^{(1)}}{\partial
u}\right\rvert_{\bm{U}^{m,i}}+\left(\bm{I}_{ij}^{m}\tilde{\bm{g}}_{v}^{m,j}\right)\left.\frac{\partial\tau_{ij}^{(2)}}{\partial
u}\right\rvert_{\bm{U}^{m,i}}\right)\left.\psi_{j}\right\rvert_{\bm{X}^{m,i}},$
(3.32)
$\displaystyle\tilde{\bm{g}}_{v}^{m,i}=\tilde{\bm{g}}_{v}^{m,i}\left.\psi_{i}\right\rvert_{\bm{X}^{m,i}}+\sum_{j\neq
i,1\leq j\leq
n_{p}}\left(\left(\bm{I}_{ij}^{m}\tilde{\bm{g}}_{u}^{m,j}\right)\left.\frac{\partial\tau_{ij}^{(1)}}{\partial
v}\right\rvert_{\bm{U}^{m,i}}+\left(\bm{I}_{ij}^{m}\tilde{\bm{g}}_{v}^{m,j}\right)\left.\frac{\partial\tau_{ij}^{(2)}}{\partial
v}\right\rvert_{\bm{U}^{m,i}}\right)\left.\psi_{j}\right\rvert_{\bm{X}^{m,i}},$
(3.33)

for $i=1,\ldots,n_{p}$, where $\bm{I}_{ij}^{m}$ is the cubic splines
interpolation operator. It returns approximated function values on
$\bm{U}^{m,i}$ grid from the known values on the grid $\bm{U}^{m,j}$.

The three steps, namely the patch extension, the finite differences and and
the blending of derivatives described above constitute our numerical scheme to
calculate the interfacial forces using the formulas given in Section A.2.

_Work complexity:_ The finite difference operators $\bm{D}_{u}^{m,i}$ and
$\bm{D}_{v}^{m,i}$ and the inter-patch interpolation operators
$\bm{I}_{ij}^{m},i,j\in\\{1,\ldots,n_{p}\\}$, can all be precomputed and
stored for application as sparse matrices. The work complexity of patch
extension is of the order of number of additional points in the extended grid,
_i.e.,_ $O(n_{p}(m+5)^{2}-n_{p}(m-1)^{2})\equiv O(n_{p}m)$. The work required
for applying finite difference operators is $O(n_{p}m^{2})$. Finally, the work
required for blending derivatives across patches is $O(n_{p}m^{2})$. Thus, the
total work complexity for calculating derivatives is $O(n_{p}m^{2})$. In our
implementation $n_{p}=6$, so the work complexity is $O(6m^{2})$.

_Accuracy:_ While the finite difference operators we use are fifth order
accurate, our numerical scheme is expected to be fourth order accurate because
we use cubic splines for patch extension and blending. This is confirmed in
our numerical experiments as tabulated in Table 2(b) and Table 2(b) and
discussed in later sections.

_Remark 3:_ We mention that the main purpose of blending is to get consistent
unique derivative values by taking weighted average across different patches
since a point $\bm{x}\in\gamma$ can belong to multiple patches. We note here
that the blending also improves the accuracy of the derivatives as evidenced
in the results tabulated in Table 7 and discussed in Section A.4.

###  3.6 Time stepping and the overall algorithm

Here, we describe the time stepping scheme we use to solve the boundary
integral formulation given by Equations 2.9, 2.10 and 2.11. The equations can
be reformulated into an initial value problem as

$\displaystyle\frac{\partial\bm{x}}{\partial t}$
$\displaystyle=\bm{u}_{\infty}(\bm{x})+\mathcal{S}_{\gamma}[\bm{f}](\bm{x}),$
(3.34) $\displaystyle=\mathcal{L}(t,\bm{x}),\ \forall\bm{x}\in\gamma,$ (3.35)

where the initial position and initial deformation gradient tensor is known.
We use an explicit Runge-Kutta-Fehlberg 4(5) [14, 15] adaptive time stepping
scheme to simulate the time evolution of capsules. This scheme is fourth-order
accurate in time. We use fixed relative tolerance $\epsilon_{tol}=10^{-6}$ for
adaptive time stepping in our simulations. Below we give the overall algorithm
for the evaluation of the right hand side of Equation 3.35:

  1. 1.

Given the initial capsule surface, we first initialize the atlas, the
transition maps and its derivatives (for Equation 3.26 and Equation 3.27). We
precompute
$\bm{\Psi}^{m},\bm{\Psi}_{u}^{m},\bm{I}_{u}^{m},\bm{I}_{d}^{m},\bm{I}_{i,j}^{m},\bm{I}_{ij}^{m,ext},\\\
\bm{D}_{u}^{m,i}$, and $\bm{D}_{v}^{m,i}$ for the $m_{\mathrm{th}}$-order
grid. We also precompute the reference tangents and normals at the
discretization points using the derivative scheme in Section 3.5. Let us
denote these reference tangents by $\bm{X}_{\partial u}^{r,m},\bm{X}_{\partial
v}^{r,m}$ and the normals by $\bm{N}^{r,m}$ respectively.

  2. 2.

At a time instant $t$, we have the discretization points $\bm{X}^{m}$ on the
capsule surface . We compute the tangents and the normals at all
discretization points in the current configuration using the scheme in Section
3.5. Let us denote these by $\bm{X}_{\partial u}^{m},\bm{X}_{\partial v}^{m}$
and $\bm{N}^{m}$.

  3. 3.

Given the reference and the current tangents along with the normals, we
compute the coefficients of the first fundamental form, the surface area
element $\bm{W}^{m}$ and the shear stress tensor $\bm{\Lambda}^{m}$.

  4. 4.

We compute the interfacial force
$\bm{F}^{m}=\nabla_{\gamma}\cdot\bm{\Lambda}^{m}$ using the derivative scheme
in Section 3.5.

  5. 5.

We compute the upsampled quantities
$\bm{X}^{m}_{u}=\bm{I}^{m}_{u}\bm{X}^{m},\bm{F}^{m}_{u}=\bm{I}_{u}^{m}\bm{F}^{m},\bm{W}^{m}_{u}=\bm{I}_{u}^{m}\bm{W}^{m}$
and the regularization parameters $\bm{\delta}^{m}_{u}$.

  6. 6.

We compute the Stokes single layer potential on the upsampled grid, denoted by
$\bm{S}^{m}_{u}$, using $\bm{X}^{m}_{u},\bm{F}^{m}_{u},\bm{\Psi}_{u}^{m}$ and
$\bm{W}^{m}_{u}$.

  7. 7.

We downsample the Stokes potential $\bm{S}^{m}=\bm{I}_{d}^{m}\bm{S}^{m}_{u}$
and add $\bm{u}_{\infty}(\bm{X}^{m})$ to get the discretized version of the
RHS of Equation 3.35 at time instant $t$ and position $\bm{X}^{m}$.

_Accuracy and work complexity of the overall algorithm:_ Our quadrature
scheme, differentiation scheme and the time stepping scheme are all fourth
order convergent. Hence, our overall scheme is fourth order convergent. We
observe this convergence numerically in the convergence results in Table 3(b).
The work complexity of our quadrature scheme and the differentiation scheme is
$O(m^{4})$ and $O(m^{2})$ respectively for $m_{\mathrm{th}}$-order grid. Thus,
the overall work complexity of our algorithm for a single time step is
$O(m^{4})$. ,

_Computation of $||\nabla_{\mathbb{S}^{2}}\phi||_{\infty}$_: In the results
section, we monitor the norm of the surface gradient of
$\phi:\mathbb{S}^{2}\longrightarrow\gamma$ with respect to the unit sphere
$\mathbb{S}^{2}$ which measures the smoothness of the capsule surface during
the simulation. It serves as a proxy for measuring the stability of capsule
dynamics during the simulation. This surface gradient
$\nabla_{\mathbb{S}^{2}}\phi$ is calculated using the surface gradient formula
given in the Section A.2 by using the first fundamental form coefficients
$E,F,G$ for the unit sphere and the surface area element $W$ of the unit
sphere $\mathbb{S}^{2}$. We use the derivative scheme in Section 3.5 to
calculate the first fundamental form using the stored discrete points
$\bm{X}^{0,m}$ on the sphere and noting that $\phi(\bm{X}^{0,m})=\bm{X}^{m}$.

##  4 Results

Now we discuss the various results to verify the accuracy and convergence of
our numerical schemes. We also present the results of simulation of extensible
capsule under shear flow and the Poiseuille flow. We validate our simulations
with the existing results in the literature. We also use the spherical
harmonics based numerical scheme used in [34] for quantitative comparisons
with our simulations. Below we give a quick summary of the spherical harmonics
based scheme used in [34].

_Spherical harmonics:_ Spherical harmonics provide an orthonormal basis [26]
for the square-integrable functions defined on the sphere $\mathbb{S}^{2}$.
Hence, they provide a spectral representation of the surfaces and have been
used for simulating Stokesian particulate flows in [34]. Below we summarize
the number of discretization points used and the work complexity for the
differentiation and singular integration in this scheme.

  1. 1.

_Discretization:_ We denote by $p$ the degree of spherical harmonics used to
represent the surface and surface fields. For degree $p$, the scheme uses
$2p(p+1)$ number of discretization points.

  2. 2.

_Singular quadrature:_ The scheme uses Graham-Sloan quadrature [17] to compute
the Stokes layer potential which is $O(p^{5})$ in work complexity.

  3. 3.

_Differentiation:_ For surface derivatives, the scheme uses spherical
harmonics based spectral differentiation which has work complexity $O(p^{3})$.
Thus, the overall work complexity of the algorithm is $O(p^{5})$.

###  4.1 Integration results

Here, we discuss the numerical accuracy and convergence of our integration
schemes discussed in Section 3.3 and Section 3.4. We present the relative
error in computing Stokes single layer potential with and without upsampling
on two different surfaces. First, we take an ellipsoid surface given by

$\displaystyle\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1,$
(4.1)

where $a=0.6,b=1,c=1$. The relative error in a quantity $q$, denoted by
$\epsilon_{q}$, is defined as
$\epsilon_{q}:=\frac{||q-q_{\mathrm{ref}}||_{\infty}}{||q_{\mathrm{ref}}||_{\infty}}$,
where $q_{\mathrm{ref}}$ is the reference solution. We use spherical harmonic
solutions with degree $p=64$ as the high accuracy reference solutions for the
surface derivatives [34]. We report the relative errors in computing single
layer Stokes potential with and without upsampling, denoted by
$\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ and $\epsilon_{\mathcal{S}[\bm{f}]}$
respectively, in Table 1(b).

As a second example, we take a 4-bump shape which can be written in standard
spherical parameters as

$\displaystyle\bm{X}(u,v)=\rho(u,v)\begin{pmatrix}\sin{u}\cos{v}\\\
\sin{u}\sin{v}\\\ \cos{u}\\\ \end{pmatrix},\ \forall u\in[0,\pi]\times
v\in[0,2\pi)$ (4.2)

where $\rho(u,v)=1+e^{-3Re(Y_{lm}(u,v))},\textit{ with }l=3,m=2$. This 4-bump
shape is shown in Figure 7. We report the relative errors in Stokes single
layer potential on this surface in Table 1(b).

The upsampling improves the accuracy of the integration scheme and helps us to
long time-horizon simulations that we discuss in the later sections. As
evident from the tabulated results, upsampling for singular quadrature is
required to get similar digits of accuracy as the derivative accuracy 111The
digits of accuracy for mean curvature $H$ and Gaussian curvature $K$ are lower
than the upsampled quadrature but since they are not required in shear force
calculation, we don’t need upsampling for derivatives. However, since the
bending forces [34] require curvature calculations, upsampling for derivatives
could be desirable if bending force is to be included.  ( Table 2(b) and Table
2(b)). The reference solution is computed using the Graham-Sloan quadrature
[17, 34] for spherical harmonics order $p=64$. We also present the relative
errors in computing area $A$ and volume $V$222Volume $V$ of surface $\gamma$
is given by $V=\int_{enc(\gamma)}dV$, where $enc(\gamma)$ refers to the volume
enclosed by $\gamma$. Using divergence theorem, we can write it as a smooth
surface integral as $V=\int_{\gamma}\bm{M}(x,y,z)\cdot\bm{n}d\gamma$, where
$\bm{M}(x,y,z)=(x,0,0)$.  to demonstrate the convergence and accuracy of our
integration scheme for smooth functions. For further verification and to study
the change in relative errors with the reduced volume $\nu$, we provide the
error in our upsampled quadrature scheme in Section A.6. In Section A.5, we
also tabulate the relative errors demonstrating the convergence of our
derivative and integration schemes for the complex shapes obtained in actual
shear and Poiseuille flow simulations discussed in the later sections. This
further verifies the correctness of our code.

Table 1: Relative error in the computation of Stokes single layer potential
(with density $\bm{f}=(x^{2},y^{2},z^{2})$ on (a) an ellipsoid and (b) the
4-bump shape (as shown in Figure 7). $\epsilon_{\mathcal{S}[\bm{f}]}$ is
relative error with no upsampling and $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ is
the relative error with 4$\times$ upsampling. $\epsilon_{A}$ and
$\epsilon_{V}$ are the relative errors in computing area and volume
respectively. $N$ is the number of discretization points. Reference solution
is computed using $p=64$ spherical harmonics (see [34]).

$m$ | $N$ | $\epsilon_{\mathcal{S}[\bm{f}]}$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{A}$ | $\epsilon_{V}$
---|---|---|---|---|---
$8$ | $294$ | $6\text{\times}{10}^{-2}$ | $7\text{\times}{10}^{-3}$ | $5\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-3}$
$16$ | $1350$ | $8\text{\times}{10}^{-3}$ | $4\text{\times}{10}^{-4}$ | $5\text{\times}{10}^{-4}$ | $4\text{\times}{10}^{-4}$
$32$ | $5766$ | $4\text{\times}{10}^{-4}$ | $2\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-5}$
$64$ | $23\,814$ | $1\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-6}$ | $1\text{\times}{10}^{-6}$ | $1\text{\times}{10}^{-6}$

(a) Relative error in Stokes single layer potential on ellipsoid of reduced
volume $\nu=0.9$, given
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.6,b=1,c=1$.

$m$ | $N$ | $\epsilon_{\mathcal{S}[\bm{f}]}$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{A}$ | $\epsilon_{V}$
---|---|---|---|---|---
$8$ | $294$ | $5\text{\times}{10}^{-1}$ | $4\text{\times}{10}^{-1}$ | $2\text{\times}{10}^{-1}$ | $2\text{\times}{10}^{-1}$
$16$ | $1350$ | $2\text{\times}{10}^{-1}$ | $9\text{\times}{10}^{-2}$ | $3\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$
$32$ | $5766$ | $3\text{\times}{10}^{-2}$ | $8\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-3}$
$64$ | $23\,814$ | $3\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-4}$ | $2\text{\times}{10}^{-4}$ | $2\text{\times}{10}^{-4}$

(b) Relative error in Stokes single layer potential for the 4-bump shape.

###  4.2 Derivative results

Now we discuss the results showing the accuracy and convergence of our
numerical differentiation scheme. We report the relative errors in derivative
calculations on the ellipsoid and the 4-bump shape in Table 2(b) and Table
2(b) respectively. The derivatives converges as the the grid length $h_{m}$
decreases. Empirically, the convergence is of the fourth order which is what
we predicted because of the cubic spline interpolation used in patch extension
operators and blending. We also provide verification of derivatives for
ellipsoids of different reduced volumes and tabulate the results in Section
A.6.

Figure 7: The 4-bump shape with its first three patches shaded with dark mesh.
Table 2: Relative error in the derivative scheme for (a) an ellipsoid and (b)
the 4-bump shape (shown in Figure 7). We report relative error in surface
normals $\bm{n}$, mean curvature $H$, Gaussian curvature $K$ and surface
divergence $div_{\gamma}$. $N$ is the number of discretization points.
Spherical harmonics results for $p=64$, using the algorithms in [34], are used
as true values and the error is computed relative to those values. Surface
divergence is computed for the smooth function
$\bm{g}(x,y,z)=(x^{2},y^{2},z^{2})$ on the surface.

m | $N$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---
$8$ | $294$ | $5\text{\times}{10}^{-3}$ | $4\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-3}$ | $4\text{\times}{10}^{-2}$
$16$ | $1350$ | $2\text{\times}{10}^{-4}$ | $5\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-4}$ | $3\text{\times}{10}^{-3}$
$32$ | $5766$ | $1\text{\times}{10}^{-5}$ | $3\text{\times}{10}^{-5}$ | $6\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-4}$
$64$ | $23\,814$ | $7\text{\times}{10}^{-7}$ | $2\text{\times}{10}^{-6}$ | $3\text{\times}{10}^{-6}$ | $1\text{\times}{10}^{-5}$

(a) Relative error in derivatives for an ellipsoid of reduced volume
$\nu=0.9$, given by
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.6,b=1,c=1$.

m | $N$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---
$8$ | $294$ | $2\text{\times}{10}^{-1}$ | $3\text{\times}{10}^{-1}$ | $8\text{\times}{10}^{-1}$ | $6\text{\times}{10}^{-1}$
$16$ | $1350$ | $6\text{\times}{10}^{-2}$ | $1\text{\times}{10}^{-1}$ | $1\text{\times}{10}^{-1}$ | $3\text{\times}{10}^{-1}$
$32$ | $5766$ | $6\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-2}$ | $6\text{\times}{10}^{-2}$
$64$ | $23\,814$ | $5\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-3}$ | $5\text{\times}{10}^{-3}$

(b) Relative error in derivatives for the 4-bump shape.

###  4.3 Accuracy and convergence of the full numerical scheme

Now, we discuss the accuracy and convergence of our full numerical solver for
the extensible capsule simulation. To this end, we take an ellipsoidal capsule
with an initial surface configuration given by
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.9,b=1,c=1$. The initial shape is taken to be the stress-free reference
configuration. We do the simulations under two different imposed background
flows, _i.e,_ shear flow and Poiseuille flow. The flows are described
mathematically as

$\displaystyle\bm{u}_{\infty}(x,y,z)$
$\displaystyle=\dot{\gamma}(y,0,0),\textit{ (shear flow) },$ (4.3)
$\displaystyle\bm{u}_{\infty}(x,y,z)$
$\displaystyle=(\alpha(R_{0}^{2}-y^{2}-z^{2}),0,0)\textit{ (Poiseuille flow)
},$ (4.4)

where $\dot{\gamma}$ is the shear rate of the shear flow, $\alpha$ controls
the curvature of the Poiseuille flow and $R_{0}$ is the radius of the circular
cross-section of the Poiseuille flow. We set $\dot{\gamma}=1$, $\alpha=1$ and
$R_{0}=5$ for these simulations. We use the differentiation, integration and
the time stepping schemes discussed above to simulate these setups for a fixed
time horizon $[0,T]$ using our solver. In Table 3(b), we tabulate the relative
errors in area $A$, volume $V$ and moment of inertia tensor $J$ of the final
capsule shape for different values of the discretization order $m$. We observe
fourth order convergence in the relative errors. The parachute shapes at time
$T=0.5$ for the Posieuille flow simulations for different levels of
discretization are given in Figure 8.

Table 3: Self-convergence results for an extensible capsule simulation under
background (a) shear flow (given by $\bm{u}_{\infty}(x,y,z)=(y,0,0)$) and (b)
Poiseuille flow (given by $\bm{u}_{\infty}(x,y,z)=((25-y^{2}-z^{2}),0,0)$).
The initial shape is an ellipsoid given by
$x^{2}/a^{2}+x^{2}/b^{2}+x^{2}/c^{2}=1,$ where $a=0.9,b=1.0,c=1.0$. We set the
shear modulus to be $E_{s}=2$ and the dilatation modulus to be $E_{D}=20$. We
simulate the capsule for a time horizon $[0,T]$ and report the relative errors
in area $A$, volume $V$ and moment of inertia tensor $J$ of the capsule shape
at time $T=0.5$. The reference solution for the relative error is the
numerical solution computed using the $m=48$ grid. We use Runge-Kutta-Fehlberg
time stepping scheme ( see Section 3.6) with fixed relative tolerance of
$\epsilon_{tol}=10^{-6}$.

$m$ | $\epsilon_{A}$ | $\epsilon_{V}$ | $\epsilon_{J}$
---|---|---|---
$8$ | $2\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$ | $3\text{\times}{10}^{-2}$
$16$ | $2\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-3}$
$32$ | $1\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-4}$

(a) Relative errors in the capsule dynamics driven by a shear flow.

$m$ | $\epsilon_{A}$ | $\epsilon_{V}$ | $\epsilon_{J}$
---|---|---|---
$8$ | $2\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-2}$
$16$ | $3\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-3}$
$32$ | $2\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-4}$

(b) Relative errors in the capsule dynamics driven by a background Poiseuille
flow.

(a)

(b)

(c)

(d)

Figure 8: The parachute shape under background Poiseuille flow (see Equation
4.4) obtained after time $T=0.5$ for different levels of discretization (a)
$m=8$, (b) $m=16$, (c) $m=32$ and (d) $m=48$ for $E_{s}=2,E_{D}=20$. The
initial shape is same as the stress-free reference state and is given by
$x^{2}/a^{2}+x^{2}/b^{2}+x^{2}/c^{2}=1,$ where $a=0.9,b=1.0,c=1.0$.

###  4.4 Relaxation of capsule

Now, we simulate the relaxation of a capsule to its stress-free reference
state to validate our code. We take a capsule with an initial shape of a unit
sphere. The initial shape is also the stress-free reference configuration of
the capsule. We simulate the dynamics of the capsule under a background linear
shear flow with shear rate $\dot{\gamma}=1$ for a fixed time horizon
$[0,T_{1}]$ followed by a zero background flow velocity for $[T_{1},T]$. We
report the resulting shapes at different times in Figure 9 . We observe that
once the background flow diminishes the capsule returns to the stress-free
state, _i.e.,_ a unit sphere. Using the moment of inertia tensor $J$ and
volume $V$ of the capsule, we compute the instantaneous Taylor asphericity
parameter [21] $D_{a}$ of the capsule shape defined as

$\displaystyle D_{a}=\frac{L-S}{L+S},$ (4.5) $\displaystyle\textit{ where
}S=\sqrt{\frac{J_{xx}+J_{yy}-\sqrt{(J_{xx}-J_{yy})^{2}+4J_{xy}^{2}}}{2V}}$
$\displaystyle,\textit{
}L=\sqrt{\frac{J_{xx}+J_{yy}+\sqrt{(J_{xx}-J_{yy})^{2}+4J_{xy}^{2}}}{2V}}.$
(4.6)

For a sphere, $D_{a}=0$. We plot the Taylor asphericity as a function of time
to monitor the relaxation back to the reference unit sphere shape and plot it
in Figure 9(d). As expected, we observe that the capsule relaxes back to a
unit sphere when the background flow is removed as $D_{a}$ drops back to zero.
Our results agree quantitatively with the results obtained using spherical
harmonics based scheme used in [34].

(a) $t=0$

(b) $t=T_{1}$

(c) $t=T$

(d)

Figure 9:  Shapes obtained at different times during simulation of the
relaxation of the capsule with initial spherical shape as the stress-free
reference state. We impose background shear flow
$\bm{u}_{\infty}(x,y,z)=(y,0,0)$ for time horizon $[0,T_{1}]$. This is
followed by zero background velocity from $[T_{1},T]$. The capsule relaxes
back to the stress-free reference shape as expected. We take
$m=16,T_{1}=1,T=17.5,E_{s}=2$ and $E_{D}=20$. (a) Shape at time $t=0$. (b)
Shape at time $t=T_{1}$. (c) Shape at time $t=T$. (d) The plot of Taylor
asphericity $D_{a}$ vs time $t$, where $D_{a}$ increases till $t=1$ under
background shear flow and then drops to zero during the relaxation phase when
the background flow is zero. Blue solid lines is the plot using our scheme and
red dotted lines denote is the plot using spherical harmonics based scheme
(with spherical harmonic degree $p=32$).

###  4.5 Capsule in shear flow

Now, we present results for steady shapes of extensible capsule under
background shear flow (see Equation 4.3). We start with a stress-free
spherical capsule. Under linear shear flow, the capsule is known to take a
terminal nearly-ellipsoidal shape and exhibit a stable tank treading motion
[11, 21]. The terminal shape and inclination angle of its major axes with the
flow direction depends on three key parameters, namely, the shear rate
$\dot{\gamma}$ of the flow, the shear modulus $E_{s}$ and the dilatation
modulus $E_{D}$ of the membrane. We simulate these dynamics for different
values of these parameters and plot the terminal inclination angles $\theta$
(with respect to the flow direction) in Figure 10(a). We compare our results
with the numerical results from [11] and obtain good quantitative agreement.
We also plot the evolution of Taylor asphericity $D_{a}$ with time in Figure
10(b) and observe good agreement with the numerical results obtained using the
spherical harmonics based scheme mentioned in [34]. Figures 11, 12 and 13
shows the terminal shapes of different reduced volumes obtained in shear flows
using our code. Our code is able to resolve shapes for different reduced
volumes obtained for a wide range of ratios of shear modulus $E_{s}$ and
dilatation modulus $E_{D}$.

To further demonstrate the effectiveness of four times upsampling, we plot the
norm of the gradient of mapping $\phi$ from a unit sphere to the capsule
surface $\gamma$ with time in shear flow simulations for different grid orders
$m$ in Figure 14. The gradient norm serves as a proxy for the stability of the
capsule shape. As evident from the plots, we need high numerical accuracy to
do long time horizon simulations. We observe that grid order $m=16$ with four
times upsampling is sufficient to do shear flow simulations. Lower upsampling
factors or lower grid order $m$ results in unstable gradients which blow up
over long time scales as shown in the Figure 14.

(a)

(b)

Figure 10:  Validation of shear flow simulation results computed using $m=16$.
(a) Plot of terminal inclination angle $\theta$ for different shear rates and
membrane mechanical properties ($E_{s}$ and $E_{D}$). We compare our results
with results in [11] and obtain good quantitative agreement. (b) We plot the
evolution of Taylor asphericity $D_{a}$ of the the capsule vs time and compare
our results (solid lines) with the spherical harmonics ($p=32$) based scheme
in [34] (dashed lines). We observe good quantitative agreement. For
$E_{D}=200$, the spherical harmonics code is unable to resolve the shapes and
the surface fields and therefore, we do not provide plots for those.

(a) $t=0$

(b) $t=0.25$

(c) $t=7.5$

Figure 11: Snapshots of capsule shape at different time instants for shear
flow simulation with $\dot{\gamma}=1,E_{s}=2,E_{D}=200$. (a) $t=0$, (b)
$t=0.25$ and (c) $t=7.5$. The shape shown in (c) is the terminal shape.

(a) $t=0.5$

(b) $t=2$

(c) $t=15$

Figure 12: Snapshots of capsule shape at different time instants for shear
flow simulation with $\dot{\gamma}=0.5,E_{s}=2,E_{D}=1$. (a) $t=0.5$, (b)
$t=2$ and (c) $t=15$. The shape shown in (c) is the terminal shape.

(a) $t=0.0417$

(b) $t=0.25$

(c) $t=1.25$

Figure 13: Snapshots of capsule shape at different time instants for shear
flow simulation with $\dot{\gamma}=6,E_{s}=2,E_{D}=1$. (a) $t=0.0417$, (b)
$t=0.25$ and (c) $t=1.25$. The shape shown in (c) is the terminal shape.

(a) Terminal reduced volume $\nu=0.85$.

(b) Terminal reduced volume $\nu=0.65$.

Figure 14:  Plots of the $||\nabla_{\mathbb{S}^{2}}\phi||_{\infty}$ vs time
$t$ in shear flow simulation where $\phi:\mathbb{S}^{2}\longrightarrow\gamma$
is the mapping from the unit sphere to the capsule $\gamma$. Initial shape is
stress-free unit sphere with (a) $\dot{\gamma}=1$, $E_{s}=2$ and $E_{D}=20$
and (b) $\dot{\gamma}=1.5$, $E_{s}=2$ and $E_{D}=1$. The plots show the
effectiveness of four times upsampling in doing stable long time horizon
simulations while two times upsampling fails for $m=16$.

###  4.6 Capsule in Poiseuille flow

In Figure 16, we present the terminal shapes for a capsule under Poiseuille
flow (see Equation 4.4) for different membrane elasticity parameters leading
to shapes of reduced volume as low as $\nu=0.6$. As for the shear flow
simulations, we plot the norm of the gradient of mapping $\phi$ from a unit
sphere to the capsule surface $\gamma$ with time for the Poiseuille flow
simulations for different grid orders $m$ in Figure 15. We observe that grid
order $m=32$ with four times upsampling is sufficient to do Poiseuille flow
simulations. Lower upsampling factors or lower grid order $m$ results in
unstable gradients which blow up over long time scales as shown in the Figure
15.

(a) Terminal reduced volume $\nu=0.96$.

(b) Terminal reduced volume $\nu=0.86$.

Figure 15:  Plots of the $||\nabla_{\mathbb{S}^{2}}\phi||_{\infty}$ vs time
$t$ in Poiseuille flow simulation where
$\phi:\mathbb{S}^{2}\longrightarrow\gamma$ is the mapping from the unit sphere
to the capsule $\gamma$. Initial shape is stress-free unit sphere with (a)
$\alpha=1.5$, $R_{0}=5$, $E_{s}=2$ and $E_{D}=200$ and (b) $\alpha=1.5$,
$R_{0}=5$, $E_{s}=2$ and $E_{D}=20$. Plots show that $m=32$ with four times
upsampling gives stable simulations for Poiseuille flow while lower values of
$m$ are not enough to resolves the shapes in Poiseuille flow.

(a) $\nu=0.96$

(b) $\nu=0.85$

(c) $\nu=0.60$

Figure 16:  Terminal parachute shapes of varying reduced volumes $\nu$ under
Poiseuille flow. (a) Terminal shape of reduced volume $\nu=0.96$ obtained for
Poiseuille flow simulation with $\alpha=1.5,E_{s}=2,E_{D}=200$. (b) Terminal
shape of reduced volume $\nu=0.85$ obtained for Poiseuille flow simulation
with $\alpha=1.5,E_{s}=2,E_{D}=20$. (c) Terminal shape of reduced volume
$\nu=0.6$ obtained for shear flow simulation with
$\alpha=1.5,E_{s}=2,E_{D}=1$. We use grid order $m=32$ and the Poiseuille flow
cross-section radius $R_{0}=5$ in these simulations .

###  4.7 Timing results

We compute the GPU wall clock times required for our scheme and give its
breakdown over different stages of our scheme. To this end, we take an
initially spherical capsule in stress free state (with $E_{s}=2,E_{D}=20$) and
simulate its dynamics under shear flow with shear rate $\dot{\gamma}=1$ till
time $T=0.1$. We do these simulations for different grid order $m$ and plot
the total time and its breakdown for different stages in Figure 17(b). We
observe that for $m\geq 16$, majority of the time is required to compute the
singular quadrature. In fact, for $m\geq 32$ virtually all of the time is
consumed in computing the quadrature. The wall clock times required to compute
the singular quadrature once is also plotted separately in Figure 17(a). We
also provide a comparison of the wall clock time per time step taken by our
GPU accelerated scheme with the spherical harmonics CPU code [34] in Section
A.8.

(a)

(b)

Figure 17:  (a) GPU wall clock times for computing singular layer quadrature
with four times upsampling for different grid orders $m=8,16,32,48$. (b)
Breakdown of wall clock times for shear flow simulation (with initial shape as
unit sphere over time horizon $T=0.1$) over different stages (computation of
singular quadrature, computation of interfacial forces and the rest of the
algorithm) for different grid orders $m$.

##  5 Fast multipole method based acceleration

Our quadrature does not involve a product quadrature and is amenable to
acceleration via fast multipole method (FMM) [33]. A full FMM acceleration
reduces the time complexity of our scheme to $O(N)$ allowing us to perform
high resolution simulations. We do not implement a multilevel FMM for the
purpose of this work and focus on showcasing the correctness and advantages of
our numerical scheme using a simple single level FMM acceleration for our
scheme. Based on the kernel independend FMM of [40], we describe our schme
below.

  1. 1.

We first use a $k$-means clustering algorithm to cluster together the
discretization points on the capsule surface for a reasonable value of $k$. We
will have $k$ clusters. We found that $k=100$ gives us relative accuracy of
accuracy of $10^{-5}$, which is sufficient for the simulations in this paper.

  2. 2.

Each cluster is imagined to be enclosed by a surrounding cubical box, called
an equivalent surface, with equivalent sources placed on a regular grid on the
boundary of this cubical box. Thus, we have $k$ equivalent surfaces and let
each equivalent surface have $N_{\mathrm{eq}}$ number of equivalent sources.
We also consider a slightly bigger cube surrounding this enclosing cube which
will be used as a check surface. The density on the equivalent sources lying
on the equivalent surface is estimated by equating the potential on the check
surface due to the actual point sources (_i.e.,_ discretization points) on the
part of the capsule surface lying inside that box or equivalent surface.

  3. 3.

We apply direct all pairs single layer potential calculation for the
neighboring boxes while the equivalent sources are used to calculate the
single layer potential for the discretization points in the far field (_i.e.,_
the discretization points inside non-neighboring boxes).

We present convergence and speedup results for calculating single layer on an
ellipsoid of reduced volume $\nu=0.9$ in Table 4(a). While we do not see
speedup for low values of $m$ because of FMM overhead computations, we see
about two times speedup for $m=96$ using our implementation of the single
level FMM on GPU.

Table 4: Relative error and speedup for FMM accelerated simulation (with 4x
upsampling) of a capsule suspended in Poiseuille flow as in Figure 16(c). All
pairs direct calculation for $m=32$ (with 4x upsampling) is used as the
reference solution and the relative error is computed in the moment of inertia
tensor of the capsule at the end of simulation. Time reported is wall clock
time in seconds taken by one time step of the simulation.
$t_{\mathrm{direct}}$ is direct calculation time and $t_{\mathrm{FMM}}$ is
time taken with the FMM acceleration. Number of equivalent surfaces or boxes
is fixed at 100. We vary the number of equivalent sources for an equivalent
surface, denoted by $N_{eq}$, with $m$.

m | $N$ | $N_{\mathrm{eq}}$ | $t_{\mathrm{direct}}$ | $t_{\mathrm{FMM}}$ | $\mathrm{Speedup}$ | $\epsilon_{\mathrm{FMM}}$
---|---|---|---|---|---|---
$32$ | $5766$ | $96$ | $10.78$ | $10.66$ | $1.02$ | $6\text{\times}{10}^{-3}$
$64$ | $23\,814$ | $128$ | $163.21$ | $104.20$ | $1.57$ | $4\text{\times}{10}^{-5}$
$96$ | $54\,150$ | $256$ | $905.31$ | $501.76$ | $1.81$ | $5\text{\times}{10}^{-6}$

(a)

.

##  6 Conclusions

In this work, we described a novel numerical scheme to simulate Stokesian
particulate flows. We described an overlapping patch based discretization of
the surfaces diffeomorphic to the unit sphere. Our numerical scheme uses the
regularized Stokes kernels [32] and finite differences on overset grid to
calculate the Stokes layer potential and interfacial elastic forces. We
presented a battery of results for the verification of our numerical scheme
using results from previous literature. We used our numerical scheme to
simulate an extensible capsule in suspended in Stokes flow and validated our
simulations. Our numerical scheme is a fourth order convergent scheme that is
$O(N^{2})$ in work complexity where $N$ is the number of discretization points
and can be accelerated to run with $O(N)$ work complexity using a multilevel
FMM. This is much better than the asymptotic work complexity of the spherical
harmonics based spectral scheme [34]. Our scheme also allows for independent
control over local resolution due to patch based parameterization of surface.
We used GPU acceleration to demonstrate the ability of our code to simulate
the complex shapes with high resolution.

In this paper, we implemented a single level FMM to demonstrate FMM based
speedup. A GPU based implementation of a multilevel FMM will further enable us
to do highly accurate simulations in a reasonable time and help us better
study the Stokesian particulate flows. We leave it as the subject of our
future work.

##  Appendix A Appendix

###  A.1 Transition maps

Here, we list the expressions for the transition maps $\tau_{ij}$ for
$i=1,..,6$, for the parameterization we use. The expressions are given as
follows:

$\displaystyle\tau_{12}(u,v)$
$\displaystyle=(u,\pi+v),\tau_{13}(u,v)=(u,\frac{\pi}{2}+v),\tau_{14}(u,v)=(u,v-\frac{\pi}{2}),$
(A.1) $\displaystyle\tau_{15}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\sin{v}),\cos^{-1}\left(\frac{\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.2) $\displaystyle\tau_{16}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\sin{v}),\cos^{-1}\left(\frac{\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right).$
(A.3)

The transition maps for the next three patches are given below:

$\displaystyle\tau_{21}(u,v)$
$\displaystyle=(u,\pi+v),\tau_{23}(u,v)=(u,v-\frac{\pi}{2}),\tau_{24}(u,v)=(u,v+\frac{\pi}{2}),$
(A.4) $\displaystyle\tau_{25}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\sin{v}),\cos^{-1}\left(\frac{-\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.5) $\displaystyle\tau_{26}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\sin{v}),\cos^{-1}\left(\frac{-\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.6) $\displaystyle\tau_{31}(u,v)$
$\displaystyle=(u,v-\pi/2),\tau_{32}(u,v)=(u,v+\frac{\pi}{2}),\tau_{34}(u,v)=(u,v+\pi),$
(A.7) $\displaystyle\tau_{35}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\cos{v}),\cos^{-1}\left(\frac{\sin{u}\sin{v}}{\sqrt{\sin^{2}{u}\sin^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.8) $\displaystyle\tau_{36}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\cos{v}),\cos^{-1}\left(\frac{\sin{u}\sin{v}}{\sqrt{\sin^{2}{u}\sin^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.9) $\displaystyle\tau_{41}(u,v)$
$\displaystyle=(u,v+\pi/2),\tau_{42}(u,v)=(u,v-\frac{\pi}{2}),\tau_{43}(u,v)=(u,v+\pi),$
(A.10) $\displaystyle\tau_{45}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\cos{v}),\cos^{-1}\left(\frac{-\sin{u}\sin{v}}{\sqrt{\sin^{2}{u}\sin^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.11) $\displaystyle\tau_{46}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\cos{v}),\cos^{-1}\left(\frac{-\sin{u}\sin{v}}{\sqrt{\sin^{2}{u}\sin^{2}{v}+\cos^{2}{u}}}\right)\right).$
(A.12)

The transition maps for $\mathcal{P}_{5}^{0}$ and $\mathcal{P}_{6}^{0}$ are
given as follows:

$\displaystyle\tau_{51}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\sin{v}),\cos^{-1}\left(\frac{\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.13) $\displaystyle\tau_{52}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\sin{v}),\cos^{-1}\left(\frac{-\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.14) $\displaystyle\tau_{53}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\sin{v}),\cos^{-1}\left(\frac{\cos{u}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.15) $\displaystyle\tau_{54}(u,v)$
$\displaystyle=\left(\cos^{-1}(\sin{u}\sin{v}),\cos^{-1}\left(\frac{-\cos{u}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.16) $\displaystyle\tau_{54}(u,v)$ $\displaystyle=(\pi-u,v),$ (A.17)
$\displaystyle\tau_{61}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\sin{v}),\cos^{-1}\left(\frac{\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.18) $\displaystyle\tau_{62}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\sin{v}),\cos^{-1}\left(\frac{-\sin{u}\cos{v}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.19) $\displaystyle\tau_{63}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\sin{v}),\cos^{-1}\left(\frac{\cos{u}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.20) $\displaystyle\tau_{64}(u,v)$
$\displaystyle=\left(\cos^{-1}(-\sin{u}\sin{v}),\cos^{-1}\left(\frac{-\cos{u}}{\sqrt{\sin^{2}{u}\cos^{2}{v}+\cos^{2}{u}}}\right)\right),$
(A.21) $\displaystyle\tau_{65}(u,v)$ $\displaystyle=(\pi-u,v).$ (A.22)

###  A.2 Surface derivative formulas

Here, we list the formulas for the first fundamental form coefficients
$E,F,G$, the unit normal $\bm{n}$, the second fundamental form coefficients
$L,M,N$, the mean curvature $H$, the Gaussian curvature $K$, the surface
divergence (of a vector field $\bm{g}$) and the surface gradient (of a scalar
function $g$) in terms of the local parameterization
$\bm{x}(u,v):D\subset\mathbb{R}^{2}\longrightarrow\gamma$. We use these
quantities in the computation of interfacial force and the verification of
surface derivatives.

Table 5: Formulas used for computing various surface derivative quantities.  Symbol |  Definition |  Symbol |  Definition
---|---|---|---
$E$ |  $\bm{x}_{u}\cdot\bm{x}_{u}$ |  $M$ |  $\bm{x}_{uv}\cdot\bm{n}$
$F$ |  $\bm{x}_{u}\cdot\bm{x}_{v}$ |  $N$ |  $\bm{x}_{vv}\cdot\bm{n}$
$G$ |  $\bm{x}_{v}\cdot\bm{x}_{v}$ |  $H$ |  $\frac{EN-2FM+GL}{2W^{2}}$
$W$ |  $\sqrt{EG-F^{2}}$ |  $\nabla_{\gamma}g$ |  $\frac{G\bm{x}_{u}-F\bm{x}_{v}}{W^{2}}g_{u}+\frac{E\bm{x}_{v}-F\bm{x}_{u}}{W^{2}}g_{v}$
$\bm{n}$ |  $\frac{\bm{x}_{u}\times\bm{x}_{v}}{W}$ |  $\nabla_{\gamma}\cdot\bm{g}$ |  $\frac{G\bm{g}_{u}-F\bm{g}_{v}}{W^{2}}\bm{x}_{u}+\frac{E\bm{g}_{v}-F\bm{g}_{u}}{W^{2}}\bm{x}_{v}$
$L$ |  $\bm{x}_{uu}\cdot\bm{n}$ |  $K$ |  $\frac{LN-M^{2}}{W^{2}}$

###  A.3 The choice of parameter $r_{0}$ for the partition of unity

Here, we briefly discuss our choice of parameter $r_{0}$ in the construction
of partition of unity in Section 3.1. We note here that for the specific
parameterization of the unit sphere we use (see Equation 3.8), we require
$r_{0}>\frac{3\pi}{12}$ to ensure that every point on $\mathbb{S}^{2}$ belongs
to the support of $\psi_{i}^{0}$ for at least one $i\in\\{1,\ldots,6\\}$. If
this is not the case, then $\\{\psi_{i}^{0}\\}_{i=1}^{6}$ ceases to be a
partition of unity on the unit sphere. Additionally, we need
$r_{0}<\frac{\pi}{2}$ so that the support of each $\psi_{i}^{0}$ is compactly
contained in $\mathcal{P}_{i}^{0}$. We tabulate some numerical results Table 6
for the relative errors in the surface derivatives on the unit sphere for
different values of $\pi>r_{0}>\frac{3\pi}{12}$. We find that
$r_{0}=\frac{5\pi}{12}$ gives optimal accuracy among the values of $r_{0}$ we
experimented with.

Table 6: Relative errors in surface derivatives. We compare the relative errors for different values of $r_{0}$.  $m$ | $r_{0}/\pi$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$
---|---|---|---
$8$ | $5.5\text{/}12$ | $2\text{\times}{10}^{-4}$ | $1.3\text{\times}{10}^{-3}$
$8$ | $5\text{/}12$ | $2\text{\times}{10}^{-4}$ | $1.5\text{\times}{10}^{-3}$
$8$ | $4\text{/}12$ | $4\text{\times}{10}^{-4}$ | $1.6\text{\times}{10}^{-3}$
$16$ | $5.5\text{/}12$ | $8\text{\times}{10}^{-6}$ | $1.7\text{\times}{10}^{-5}$
$16$ | $5\text{/}12$ | $7\text{\times}{10}^{-6}$ | $1.8\text{\times}{10}^{-5}$
$16$ | $4\text{/}12$ | $7\text{\times}{10}^{-6}$ | $1.9\text{\times}{10}^{-5}$
$32$ | $5.5\text{/}12$ | $4\text{\times}{10}^{-7}$ | $1.3\text{\times}{10}^{-6}$
$32$ | $5\text{/}12$ | $4\text{\times}{10}^{-7}$ | $1.7\text{\times}{10}^{-6}$
$32$ | $4\text{/}12$ | $4\text{\times}{10}^{-7}$ | $1.9\text{\times}{10}^{-6}$

###  A.4 Effect of the blending process in calculating derivatives

Here, we provide the relative errors in the computing normals $n$ and the mean
curvature $H$ for the unit sphere with and without blending process discussed
in Section 3.5. The results are tabulated in Table 7. The results show that
blending improves the accuracy of derivatives.

Table 7:  Relative error in computing surface normals $n$ and the mean curvature $H$ for unit sphere with and without the blending. $\epsilon^{\mathrm{b}}$ are the errors with blending and $\epsilon^{\mathrm{nb}}$ are the errors without the blending process in computing derivatives.  $m$ | $\epsilon_{\bm{n}}^{\mathrm{b}}$ | $\epsilon_{H}^{\mathrm{b}}$ | $\epsilon_{\bm{n}}^{\mathrm{nb}}$ | $\epsilon_{H}^{\mathrm{nb}}$
---|---|---|---|---
$8$ | $2\text{\times}{10}^{-4}$ | $1.3\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-2}$
$16$ | $8\text{\times}{10}^{-6}$ | $1.7\text{\times}{10}^{-5}$ | $9\text{\times}{10}^{-5}$ | $2\text{\times}{10}^{-3}$
$32$ | $4\text{\times}{10}^{-7}$ | $1.3\text{\times}{10}^{-6}$ | $3\text{\times}{10}^{-6}$ | $3\text{\times}{10}^{-4}$

###  A.5 Surface derivative and singular integration errors for shear flow and
Poiseuille flow terminal shapes

To further verify the convergence and accuracy of our derivative and
integration schemes, we upsample the low reduced volume shapes obtained in
Figure 13(c) and Figure 16(c) to $p=64$ spherical harmonics and use the
spherical harmonics derivatives and Graham-Sloan quadrature [34] as the
reference values to study the convergence of our schemes. We report the
relative errors for these shapes in Table 8(b).

Table 8: Relative errors in the derivative scheme and singular quadrature for
the shapes in Figure 13(c) and Figure 16(c) using our scheme (at the top for
different grid orders $m$ with four times upsampling to compute quadrature ).
The shapes are upsampled to spherical harmonics degree $p=64$ and the error is
computed relative to those values. Surface divergence is computed for the
smooth function $\bm{f}(x,y,z)=(x^{2},y^{2},z^{2})$ on the surface.

m | $N$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---|---
$8$ | $294$ | $5\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-2}$ | $5\text{\times}{10}^{-1}$ | $1$ | $9\text{\times}{10}^{-2}$
$16$ | $1350$ | $2\text{\times}{10}^{-2}$ | $6\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-1}$ | $2\text{\times}{10}^{-1}$ | $5\text{\times}{10}^{-2}$
$32$ | $5766$ | $5\text{\times}{10}^{-3}$ | $9\text{\times}{10}^{-3}$ | $8\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-3}$
$64$ | $23\,814$ | $5\text{\times}{10}^{-4}$ | $9\text{\times}{10}^{-4}$ | $8\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-4}$

(a) Relative error in singular quadrature and derivatives for the shape in
Figure 13(c) of reduced volume $\nu=0.6$.

m | $N$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---|---
$8$ | $294$ | $8\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-2}$ | $7\text{\times}{10}^{-1}$ | $8\text{\times}{10}^{-1}$ | $9\text{\times}{10}^{-2}$
$16$ | $1350$ | $2\text{\times}{10}^{-2}$ | $5\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-1}$ | $3\text{\times}{10}^{-1}$ | $6\text{\times}{10}^{-2}$
$32$ | $5766$ | $3\text{\times}{10}^{-3}$ | $8\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-2}$ | $7\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-3}$
$64$ | $23\,814$ | $6\text{\times}{10}^{-4}$ | $9\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-2}$ | $1\text{\times}{10}^{-2}$ | $8\text{\times}{10}^{-4}$

(b) Relative error in singular quadrature and derivatives for the shape Figure
16(c) of reduced volume $\nu=0.60$.

###  A.6 Numerical errors for different reduced volume $\nu$

Here, we tabulate errors in computing the singular quadrature and surface
derivatives using our numerical scheme. We report these errors for ellipsoids
of varying reduced volumes $\nu\in\\{0.92,0.78,0.5,0.3\\}$ in Table 9(b) and
Table 10(b). The tables demonstrate that the accuracy deteriorates with
decreasing reduced volumes. In general, it is harder to do long time horizon
simulations of capsules with lower reduced volumes due to this reason.

Table 9: Relative errors in the derivative scheme and singular quadrature for
ellipsoids of reduced volume (a) $\nu=0.92$ and (b) $\nu=0.78$ using our
scheme (at the top for different grid orders $m$ with four times upsampling to
compute quadrature ) and using spherical harmonics scheme [34] (at the bottom
for different spherical harmonics degrees $p$). Spherical harmonics results
for $p=64$ are used as true values and the error is computed relative to those
values. Surface divergence is computed for the smooth function
$\bm{f}(x,y,z)=(x^{2},y^{2},z^{2})$ on the surface.

m | $N$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---|---
$8$ | $294$ | $7\text{\times}{10}^{-3}$ | $5\text{\times}{10}^{-3}$ | $4\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-3}$ | $4\text{\times}{10}^{-2}$
$16$ | $1350$ | $4\text{\times}{10}^{-4}$ | $2\text{\times}{10}^{-4}$ | $5\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-3}$
$32$ | $5766$ | $2\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-5}$ | $3\text{\times}{10}^{-5}$ | $6\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-4}$
$64$ | $23\,814$ | $1\text{\times}{10}^{-6}$ | $7\text{\times}{10}^{-7}$ | $2\text{\times}{10}^{-6}$ | $3\text{\times}{10}^{-6}$ | $1\text{\times}{10}^{-5}$
p | $N$ | $\epsilon_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
$8$ | $144$ | $4\text{\times}{10}^{-5}$ | $3\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-3}$
$16$ | $544$ | $3\text{\times}{10}^{-8}$ | $1\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-5}$ | $4\text{\times}{10}^{-5}$ | $2\text{\times}{10}^{-5}$
$24$ | $1200$ | $5\text{\times}{10}^{-11}$ | $4\text{\times}{10}^{-8}$ | $6\text{\times}{10}^{-8}$ | $2\text{\times}{10}^{-7}$ | $1\text{\times}{10}^{-7}$
$32$ | $2112$ | $6\text{\times}{10}^{-13}$ | $5\text{\times}{10}^{-11}$ | $8\text{\times}{10}^{-10}$ | $4\text{\times}{10}^{-9}$ | $3\text{\times}{10}^{-9}$

(a) Relative error in singular quadrature and derivatives for an ellipsoid of
reduced volume $\nu=0.9$, given by
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.6,b=1,c=1$.

m | $N$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---|---
$8$ | $294$ | $7\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$
$16$ | $1350$ | $4\text{\times}{10}^{-4}$ | $9\text{\times}{10}^{-4}$ | $3\text{\times}{10}^{-3}$ | $5\text{\times}{10}^{-3}$ | $9\text{\times}{10}^{-4}$
$32$ | $5766$ | $2\text{\times}{10}^{-5}$ | $8\text{\times}{10}^{-5}$ | $3\text{\times}{10}^{-4}$ | $5\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-4}$
$64$ | $23\,814$ | $1\text{\times}{10}^{-6}$ | $4\text{\times}{10}^{-6}$ | $1\text{\times}{10}^{-5}$ | $3\text{\times}{10}^{-5}$ | $9\text{\times}{10}^{-6}$
p | $N$ | $\epsilon_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
$8$ | $144$ | $7\text{\times}{10}^{-4}$ | $2\text{\times}{10}^{-2}$ | $3\text{\times}{10}^{-2}$ | $6\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$
$16$ | $544$ | $2\text{\times}{10}^{-6}$ | $6\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-3}$
$24$ | $1200$ | $3\text{\times}{10}^{-8}$ | $2\text{\times}{10}^{-5}$ | $5\text{\times}{10}^{-5}$ | $1\text{\times}{10}^{-4}$ | $5\text{\times}{10}^{-5}$
$32$ | $2112$ | $5\text{\times}{10}^{-10}$ | $1\text{\times}{10}^{-7}$ | $6\text{\times}{10}^{-7}$ | $9\text{\times}{10}^{-6}$ | $4\text{\times}{10}^{-6}$

(b) Relative error in singular quadrature and derivatives for an ellipsoid of
reduced volume $\nu=0.78$, given by
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.4,b=1,c=1$.  Table 10: Relative errors in the derivative scheme and
singular quadrature for ellipsoids of reduced volume (a) $\nu=0.5$ and (b)
$\nu=0.3$ using our scheme (at the top for different grid orders $m$ with four
times upsampling to compute quadrature ) and using spherical harmonics scheme
[34] (at the bottom for different spherical harmonics degrees $p$). Spherical
harmonics results for $p=64$ are used as true values and the error is computed
relative to those values. Surface divergence is computed for the smooth
function $\bm{f}(x,y,z)=(x^{2},y^{2},z^{2})$ on the surface.

m | $N$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---|---
$8$ | $294$ | $2\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-1}$ | $2\text{\times}{10}^{-1}$ | $3\text{\times}{10}^{-1}$ | $4\text{\times}{10}^{-2}$
$16$ | $1350$ | $4\text{\times}{10}^{-3}$ | $9\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-2}$ | $6\text{\times}{10}^{-2}$ | $3\text{\times}{10}^{-3}$
$32$ | $5766$ | $7\text{\times}{10}^{-4}$ | $8\text{\times}{10}^{-4}$ | $4\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-4}$
$64$ | $23\,814$ | $6\text{\times}{10}^{-5}$ | $9\text{\times}{10}^{-5}$ | $4\text{\times}{10}^{-4}$ | $7\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-5}$
p | $N$ | $\epsilon_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
$8$ | $144$ | $9\text{\times}{10}^{-3}$ | $9\text{\times}{10}^{-2}$ | $3\text{\times}{10}^{-1}$ | $4\text{\times}{10}^{-1}$ | $2\text{\times}{10}^{-3}$
$16$ | $544$ | $2\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-2}$ | $6\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$
$24$ | $1200$ | $1\text{\times}{10}^{-5}$ | $3\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-2}$ | $7\text{\times}{10}^{-3}$
$32$ | $2112$ | $2\text{\times}{10}^{-6}$ | $1\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-3}$ | $5\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-4}$

(a) Relative error in singular quadrature and derivatives for an ellipsoid of
reduced volume $\nu=0.5$, given by
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.2,b=1,c=1$.

m | $N$ | $\epsilon^{up}_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
---|---|---|---|---|---|---
$8$ | $294$ | $4\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-1}$ | $5\text{\times}{10}^{-1}$ | $6\text{\times}{10}^{-1}$ | $4\text{\times}{10}^{-2}$
$16$ | $1350$ | $1\text{\times}{10}^{-2}$ | $9\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-1}$ | $3\text{\times}{10}^{-1}$ | $3\text{\times}{10}^{-3}$
$32$ | $5766$ | $4\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$ | $7\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-3}$
$64$ | $23\,814$ | $1\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-3}$ | $8\text{\times}{10}^{-3}$ | $9\text{\times}{10}^{-3}$ | $9\text{\times}{10}^{-4}$
p | $N$ | $\epsilon_{\mathcal{S}[\bm{f}]}$ | $\epsilon_{\bm{n}}$ | $\epsilon_{H}$ | $\epsilon_{K}$ | $\epsilon_{div_{\gamma}}$
$8$ | $144$ | $4\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-1}$ | $1$ | $1$ | $4\text{\times}{10}^{-1}$
$16$ | $544$ | $3\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-1}$ | $3\text{\times}{10}^{-1}$ | $4\text{\times}{10}^{-1}$ | $2\text{\times}{10}^{-1}$
$24$ | $1200$ | $5\text{\times}{10}^{-4}$ | $4\text{\times}{10}^{-2}$ | $1\text{\times}{10}^{-1}$ | $1\text{\times}{10}^{-1}$ | $6\text{\times}{10}^{-2}$
$32$ | $2112$ | $4\text{\times}{10}^{-5}$ | $5\text{\times}{10}^{-3}$ | $3\text{\times}{10}^{-2}$ | $4\text{\times}{10}^{-2}$ | $5\text{\times}{10}^{-3}$

(b) Relative error in singular quadrature and derivatives for an ellipsoid of
reduced volume $\nu=0.3$ given by
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.1,b=1,c=1$.

###  A.7 Sensitivity of the singular quadrature scheme to the regularization
parameter $\delta$

Here, we discuss the effect of choice of regularization parameter $\delta$ on
the singular quadrature accuracy. We experiment with different values of the
constant in $C$ in Equation 3.20. We tabulate the error in singular quadrature
for different values of $C$ in Table 11(a). We find out that $C=1$ works best
for the parameterization we have chosen. We also tabulate in the same table
the results for fixed global regularization parameter $\delta=Ch_{m}$ and show
that our patch-dependent regularization parameter gives better accuracy and
convergence compared to the fixed regularization parameter.

Table 11: Sensitivity of singular quadrature with respect to the regularization parameter $\delta$. First three columns show the relative errors for the patch-dependent $\delta=C\delta^{*}$ where $C$ is a constant. The last three columns show the relative errors for fixed global regularization parameter $\delta=Ch_{m}$.  $m$ | $\delta=0.5\delta^{*}$ | $\delta=\delta^{*}$ | $\delta=2\delta^{*}$ | $\delta=0.5h_{m}$ | $\delta=h_{m}$ | $\delta=2h_{m}$
---|---|---|---|---|---|---
$8$ | $2\text{\times}{10}^{-2}$ | $7\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-2}$ | $2\text{\times}{10}^{-2}$ | $7\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-2}$
$16$ | $6\text{\times}{10}^{-3}$ | $4\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-3}$ | $6\text{\times}{10}^{-3}$ | $5\text{\times}{10}^{-4}$ | $5\text{\times}{10}^{-3}$
$32$ | $3\text{\times}{10}^{-3}$ | $2\text{\times}{10}^{-5}$ | $6\text{\times}{10}^{-5}$ | $4\text{\times}{10}^{-3}$ | $1\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-3}$
$64$ | $2\text{\times}{10}^{-4}$ | $1\text{\times}{10}^{-6}$ | $6\text{\times}{10}^{-6}$ | $7\text{\times}{10}^{-4}$ | $6\text{\times}{10}^{-5}$ | $6\text{\times}{10}^{-4}$
(a) Relative error in the computation of Stokes single layer potential (with
density $\bm{f}=(x^{2},y^{2},z^{2})$) on the ellipsoid
$\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}+\frac{z^{2}}{c^{2}}=1$, where
$a=0.4,b=1,c=1$, for different choices of regularization parameter $\delta$.
Reference solution is computed using $p=64$ spherical harmonics.

###  A.8 Wall clock times for our GPU accelerated code _vs_ the spherical
harmonics CPU code

We report the wall clock time per time step required for our GPU-accelerated
code compared to the spherical harmonics CPU code in Table 12. We note that
the spherical harmonics discretization of degree $p$ contains $2p(p+1)$ points
compared to $6(m-1)^{2}$ for our scheme with $m_{\mathrm{th}}$-order grid.
Since, we use four times upsampling, the number of discretization points for
evaluating singular quadrature is even higher at $N_{\mathrm{up}}=6(4m-1)^{2}$
in our scheme. Even though our discretization has much higher points than the
spherical harmonics, our scheme has lower runtime for higher order grids and
hence, scales better than the spherical harmonics as the number of points
increase. Since the quadrature scheme requires most of the time in our
simulations Figure 17(b), our scheme can be further accelerated using fast
multipole methods(FMMs) [33] and can be used to do faster simulations at
higher resolutions. We leave the FMM acceleration as the subject for future
work.

Table 12: Wall clock time per time step (denoted by $T_{\mathrm{step}}$ in _seconds_) for our code with 4x upsampling (on the left) vs the spherical harmonics scheme (on the right). $N$ denotes the number of discretization points. $N_{\mathrm{up}}$ is the upsampled number of discretization points used for the singular quadrature in our scheme.  $m$ | $N$ | $N_{\mathrm{up}}$ | $T_{\mathrm{step}}$ | $p$ | $N$ | $T_{\mathrm{step}}$
---|---|---|---|---|---|---
$8$ | $294$ | $4704$ | $0.61$ | $8$ | $144$ | $0.75$
$16$ | $1350$ | $21\,600$ | $1.26$ | $16$ | $544$ | $1.92$
$32$ | $5766$ | $92\,256$ | $10.10$ | $32$ | $2112$ | $10.05$
$48$ | $13\,254$ | $381\,024$ | $46.46$ | $64$ | $8320$ | $87.86$

## References

  * [1] D. Agarwal and G. Biros.  Stable shapes of three-dimensional vesicles in unconfined and confined poiseuille flow.  Phys. Rev. Fluids, 5:013603, 2020.
  * [2] D. Agarwal and G. Biros.  Shape dynamics of a red blood cell in poiseuille flow.  Phys. Rev. Fluids, 7(9):093602, 2022.
  * [3] P. Bagchi.  Mesoscale simulation of blood flow in small vessels.  Biophs. J., 92(6):1858–1877, 2007.
  * [4] J. Beale, W. Ying, and J. Wilson.  A simple method for computing singular or nearly singular integrals on closed surfaces.  Commun. Comput. Phys., 20(3):733–753, 2016.
  * [5] T. Biben, A. Farutin, and C. Misbah.  Three-dimensional vesicles under shear flow: numerical study of dynamics and phase diagram.  Phys. Rev. E, 83:031921, 2011.
  * [6] G. Boedec, M. Leonetti, and M. Jaeger.  3d vesicle dynamics simulations with a linearly triangulated surface.  J. Comput. Phys., 230:1020–1034, 2011.
  * [7] O. Bruno and L. Kunyansky.  A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications.  Journal of Computational Physics, 169:80–110, 2001.
  * [8] D. Cordasco, A. Yazdani, and P. Bagchi.  Comparison of erythrocyte dynamics in shear flow under different stress-free configurations.  Phys. Fluids, 26:041902, 2014.
  * [9] C. de Boor.  A practical guide to splines.  Springer-Verlag, 1978.
  * [10] M. Delfour and J. Zolésio.  Shapes and geometries: metrics, analysis, differential calculus, and optimization.  2011\.
  * [11] A. Farutin, T. Biben, and C. Misbah.  3D numerical simulations of vesicle and inextensible capsule dynamics.  Journal of Computational Physics, 275:539, 2014.
  * [12] D. A. Fedosov, W. Pan, B. Caswell, G. Gompper, and G. E. Karniadakis.  Predicting human blood viscosity in silico,.  Proc. Natl. Acad. Sci., 108:11772, 2011.
  * [13] D. A. Fedosov, M. Peltomäki, and G. Gompper.  Deformation and dynamics of red blood cells in flow through cylindrical microchannels.  Soft Matter, 10:4258, 2014.
  * [14] E. Fehlberg.  Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control.  NASA Technical Report, 287, 1968.
  * [15] E. Fehlberg.  Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems.  National aeronautics and space administration, 315, 1969.
  * [16] Fornberg and Bengt.  Generation of finite difference formulas on arbitrarily spaced grids.  Mathematics of Computation, 51(184):699–706, 1988.
  * [17] I. Graham and I. Sloan.  Fully discrete spectral boundary integral methods for Helmholtz problems on smooth closed surfaces in $\mathbb{R}^{3}$.  Numerische Mathematik, 92, 2002.
  * [18] L. Greengard, M. O’Neil, M. Rachh, and F. Vico.  Fast multipole methods for the evaluation of layer potentials with locally-corrected quadratures.  Journal of Computational Physics: X, 10:100092, 2021.
  * [19] A. Guckenberger, A. Kihm, T. John, C. Wagner, and S. Gekle.  Numerical-experimental observation of shape bistability of red blood cells flowing in a microchannel.  Soft Matter, 14:2032, 2018.
  * [20] T. Kruger, M. Gross, D. Raabe, and F. Varnik.  Predicting human blood viscosity in silico,.  Soft Matter, 9:9008, 2013.
  * [21] E. Lac, D. Barthès-Biesel, N. Pelekasis, and J. Tsamopoulos.  Spherical capsules in three-dimensional unbounded Stokes flows: effect of the membrane constitutive law and onset of buckling.  J. Fluid Mech., 516:303–334, 2004.
  * [22] R. MacMeccan, J. Clausen, G. Neitzel, and C. Aidun.  Simulating deformable particle suspensions using a coupled lattice-boltzmann and finite-element method.  J. Fluid Mech., 618:13–39, 2009.
  * [23] J. Mauer, S. Mendez, L. Lanotte, F. Nicoud, M. Abkarian, G. Gompper, and D. A. Fedosov.  Flow-induced transitions of red blood cell shapes under shear.  Phys. Rev. Lett., 121:118103, 2018.
  * [24] H.-N. Nguyen and R. Cortez.  Reduction of the regularization error of the method of regularized stokeslets for a rigid object immersed in a three dimensional Stokes flow.  Commun. Comput. Phys., 15:126–152, 2014.
  * [25] L. Nylons, M. Harris, and J. Prins.  Fast N-body simulation with CUDA. GPU Gems 3.  Addison Wesley, page 677–795, 2007.
  * [26] S. Orszag.  Fourier series on spheres.  Monthly Weather Review, 102:56–75, 1974.
  * [27] C. Pozrikidis.  Boundary integral and singularity methods for linearized viscous flow.  Cambridge University Press, Cambridge, 1992.
  * [28] C. Pozrikidis.  Effect of membrane bending stiffness on the deformation of capsules in simple shear flow.  Journal of Fluid Mechanics, 440:269–291, 2001.
  * [29] C. Pozrikidis.  Interfacial dynamics for Stokes flow.  Journal of Computational Physics, 169:250–301, 2001.
  * [30] R. Skalak, A. Tozeren, R. P. Zarda, and S. Chien.  Strain energy function of red blood cell membranes.  Biophys. J., 13:245, 1973.
  * [31] S. Sukumaran and U. Seifert.  Influence of shear flow on vesicles near a wall: a numerical study.  Phys. Rev. E, 64:011916, 2001.
  * [32] S. Tlupova and J. Beale.  Regularized single and double layer integrals in 3D Stokes flow.  Journal of Computational Physics, 386:568, 2019.
  * [33] A. Tornberg and L. Greengard.  A fast multipole method for the three-dimensional Stokes equations.  J. Comput. Phys., 227:1613–1619, 2008.
  * [34] S. K. Veerapaneni, A. Rahimian, G. Biros, and D. Zorin.  A fast algorithm for simulating vesicle flows in three dimensions.  Journal of Computational Physics, 230:5610, 2011.
  * [35] V. Vitkova, M.-A. Mader, B. Polack, C. Misbah, and T. Podgorski.  Micro-Macro link in rheology of erythrocyte and vesicle suspensions?  Biophys. J., 95:L33, 2008.
  * [36] M. Wala and A. Klöckner.  A fast algorithm for quadrature by expansion in three dimensions.  Journal of Computational Physics, 388:655–689, 2019.
  * [37] F. W. Warner.  Foundations of differentiable manifolds and lie groups.  New York: Springer-Verlag, 1983.
  * [38] J. Weiner.  On a problem of Chen, Willmore, et al.  Indiana University, Mathematics Journal, 27:19–35, 1978.
  * [39] B. Wu and P.-G. Martinsson.  A unified trapezoidal quadrature method for singular and hypersingular boundary integral operators on curved surfaces.  SIAM Journal on Numerical Analysis, 61(5):2182–2208, 2023.
  * [40] L. Ying, G. Biros, and D. Zorin.  A kernel-independent adaptive fast multipole algorithm in two and three dimensions.  Journal of Computational Physics, 196:591–626, 2004.
  * [41] L. Ying, G. Biros, and D. Zorin.  A high-order 3d boundary integral equation solver for elliptic pdes in smooth domains.  Journal of Computational Physics, 219(1):247–275, 2006.
  * [42] H. Zhao, A. Isfahani, L. Olson, and J. Freund.  A spectral boundary integral method for flowing blood cells.  Journal of Computational Physics, 229:10, 2010.
  * [43] H. Zhao, A. Spann, and E. Shaqfeh.  The dynamics of a vesicle in a wall-bound shear flow.  Phys. Fluids, 23:121901, 2011.