__copyright__ = "Copyright (C) 2017 Natalie Beams"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
"""
import numpy as np
from pytential import sym
from sumpy.kernel import StokesletKernel, StressletKernel, LaplaceKernel
__doc__ = """
.. autoclass:: StokesletWrapper
.. autoclass:: StressletWrapper
.. autoclass:: StokesOperator
.. autoclass:: HsiaoKressExteriorStokesOperator
.. autoclass:: HebekerExteriorStokesOperator
"""
# {{{ StokesletWrapper
[docs]class StokesletWrapper:
"""Wrapper class for the :class:`~sumpy.kernel.StokesletKernel` kernel.
This class is meant to shield the user from the messiness of writing
out every term in the expansion of the double-indexed Stokeslet kernel
applied to the density vector. The object is created
to do some of the set-up and bookkeeping once, rather than every
time we want to create a symbolic expression based on the kernel -- say,
once when we solve for the density, and once when we want a symbolic
representation for the solution, for example.
The :meth:`apply` function returns the integral expressions needed for
the vector velocity resulting from convolution with the vector density,
and is meant to work similarly to calling
:func:`~pytential.symbolic.primitives.S` (which is
:class:`~pytential.symbolic.primitives.IntG`).
Similar functions are available for other useful things related to
the flow: :meth:`apply_pressure`, :meth:`apply_derivative` (target derivative),
:meth:`apply_stress` (applies symmetric viscous stress tensor in
the requested direction).
.. attribute:: kernel_dict
The dictionary allows us to exploit symmetry -- that
:math:`S_{01}` is identical to :math:`S_{10}` -- and avoid creating
multiple expansions for the same kernel in a different ordering.
.. automethod:: __init__
.. automethod:: apply
.. automethod:: apply_pressure
.. automethod:: apply_derivative
.. automethod:: apply_stress
"""
[docs] def __init__(self, dim=None):
self.dim = dim
if dim == 2:
self.kernel_dict = {
(2, 0): StokesletKernel(dim=2, icomp=0, jcomp=0),
(1, 1): StokesletKernel(dim=2, icomp=0, jcomp=1),
(0, 2): StokesletKernel(dim=2, icomp=1, jcomp=1)
}
elif dim == 3:
self.kernel_dict = {
(2, 0, 0): StokesletKernel(dim=3, icomp=0, jcomp=0),
(1, 1, 0): StokesletKernel(dim=3, icomp=0, jcomp=1),
(1, 0, 1): StokesletKernel(dim=3, icomp=0, jcomp=2),
(0, 2, 0): StokesletKernel(dim=3, icomp=1, jcomp=1),
(0, 1, 1): StokesletKernel(dim=3, icomp=1, jcomp=2),
(0, 0, 2): StokesletKernel(dim=3, icomp=2, jcomp=2)
}
else:
raise ValueError("unsupported dimension given to StokesletWrapper")
[docs] def apply(self, density_vec_sym, mu_sym, qbx_forced_limit):
"""Symbolic expressions for integrating Stokeslet kernel.
Returns an object array of symbolic expressions for the vector
resulting from integrating the dyadic Stokeslet kernel with
variable *density_vec_sym*.
:arg density_vec_sym: a symbolic vector variable for the density vector.
:arg mu_sym: a symbolic variable for the viscosity.
:arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on
to :class:`~pytential.symbolic.primitives.IntG`.
"""
sym_expr = np.empty((self.dim,), dtype=object)
for comp in range(self.dim):
# Start variable count for kernel with 1 for the requested result
# component
base_count = np.zeros(self.dim, dtype=np.int32)
base_count[comp] += 1
for i in range(self.dim):
var_ctr = base_count.copy()
var_ctr[i] += 1
ctr_key = tuple(var_ctr)
if i < 1:
sym_expr[comp] = sym.int_g_vec(
self.kernel_dict[ctr_key], density_vec_sym[i],
qbx_forced_limit=qbx_forced_limit, mu=mu_sym)
else:
sym_expr[comp] = sym_expr[comp] + sym.int_g_vec(
self.kernel_dict[ctr_key], density_vec_sym[i],
qbx_forced_limit=qbx_forced_limit, mu=mu_sym)
return sym_expr
[docs] def apply_pressure(self, density_vec_sym, mu_sym, qbx_forced_limit):
"""Symbolic expression for pressure field associated with the Stokeslet."""
from pytential.symbolic.mappers import DerivativeTaker
kernel = LaplaceKernel(dim=self.dim)
for i in range(self.dim):
if i < 1:
sym_expr = DerivativeTaker(i).map_int_g(
sym.int_g_vec(kernel, density_vec_sym[i],
qbx_forced_limit=qbx_forced_limit))
else:
sym_expr = sym_expr + (DerivativeTaker(i).map_int_g(
sym.int_g_vec(kernel, density_vec_sym[i],
qbx_forced_limit=qbx_forced_limit)))
return sym_expr
[docs] def apply_derivative(self, deriv_dir, density_vec_sym,
mu_sym, qbx_forced_limit):
"""Symbolic derivative of velocity from Stokeslet.
Returns an object array of symbolic expressions for the vector
resulting from integrating the *deriv_dir* target derivative of the
dyadic Stokeslet kernel with variable *density_vec_sym*.
:arg deriv_dir: integer denoting the axis direction for the derivative.
:arg density_vec_sym: a symbolic vector variable for the density vector.
:arg mu_sym: a symbolic variable for the viscosity.
:arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on
to :class:`~pytential.symbolic.primitives.IntG`.
"""
from pytential.symbolic.mappers import DerivativeTaker
sym_expr = np.empty((self.dim,), dtype=object)
for comp in range(self.dim):
# Start variable count for kernel with 1 for the requested result
# component
base_count = np.zeros(self.dim, dtype=np.int32)
base_count[comp] += 1
for i in range(self.dim):
var_ctr = base_count.copy()
var_ctr[i] += 1
ctr_key = tuple(var_ctr)
if i < 1:
sym_expr[comp] = DerivativeTaker(deriv_dir).map_int_g(
sym.int_g_vec(self.kernel_dict[ctr_key],
density_vec_sym[i],
qbx_forced_limit=qbx_forced_limit,
mu=mu_sym))
else:
sym_expr[comp] = sym_expr[comp] + DerivativeTaker(
deriv_dir).map_int_g(
sym.int_g_vec(self.kernel_dict[ctr_key],
density_vec_sym[i],
qbx_forced_limit=qbx_forced_limit,
mu=mu_sym))
return sym_expr
[docs] def apply_stress(self, density_vec_sym, dir_vec_sym,
mu_sym, qbx_forced_limit):
r"""Symbolic expression for viscous stress applied to a direction.
Returns a vector of symbolic expressions for the force resulting
from the viscous stress
.. math::
-p \delta_{ij} + \mu (\nabla_i u_j + \nabla_j u_i)
applied in the direction of *dir_vec_sym*.
Note that this computation is very similar to computing
a double-layer potential with the Stresslet kernel in
:class:`StressletWrapper`. The difference is that here the direction
vector is applied at the target points, while in the Stresslet the
direction is applied at the source points.
:arg density_vec_sym: a symbolic vector variable for the density vector.
:arg dir_vec_sym: a symbolic vector for the application direction.
:arg mu_sym: a symbolic variable for the viscosity.
:arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on
to :class:`~pytential.symbolic.primitives.IntG`.
"""
import itertools
sym_expr = np.empty((self.dim,), dtype=object)
stresslet_obj = StressletWrapper(dim=self.dim)
for comp in range(self.dim):
# Start variable count for kernel with 1 for the requested result
# component
base_count = np.zeros(self.dim, dtype=np.int32)
base_count[comp] += 1
for i, j in itertools.product(range(self.dim), range(self.dim)):
var_ctr = base_count.copy()
var_ctr[i] += 1
var_ctr[j] += 1
ctr_key = tuple(var_ctr)
if i + j < 1:
sym_expr[comp] = dir_vec_sym[i] * sym.int_g_vec(
stresslet_obj.kernel_dict[ctr_key],
density_vec_sym[j],
qbx_forced_limit=qbx_forced_limit, mu=mu_sym)
else:
sym_expr[comp] = sym_expr[comp] + dir_vec_sym[i] * sym.int_g_vec(
stresslet_obj.kernel_dict[ctr_key],
density_vec_sym[j],
qbx_forced_limit=qbx_forced_limit,
mu=mu_sym)
return sym_expr
# }}}
# {{{ StressletWrapper
[docs]class StressletWrapper:
"""Wrapper class for the :class:`~sumpy.kernel.StressletKernel` kernel.
This class is meant to shield the user from the messiness of writing
out every term in the expansion of the triple-indexed Stresslet
kernel applied to both a normal vector and the density vector.
The object is created to do some of the set-up and bookkeeping once,
rather than every time we want to create a symbolic expression based
on the kernel -- say, once when we solve for the density, and once when
we want a symbolic representation for the solution, for example.
The :meth:`apply` function returns the integral expressions needed for
convolving the kernel with a vector density, and is meant to work
similarly to :func:`~pytential.symbolic.primitives.S` (which is
:class:`~pytential.symbolic.primitives.IntG`).
Similar functions are available for other useful things related to
the flow: :meth:`apply_pressure`, :meth:`apply_derivative` (target derivative),
:meth:`apply_stress` (applies symmetric viscous stress tensor in
the requested direction).
.. attribute:: kernel_dict
The dictionary allows us to exploit symmetry -- that
:math:`T_{012}` is identical to :math:`T_{120}` -- and avoid creating
multiple expansions for the same kernel in a different ordering.
.. automethod:: __init__
.. automethod:: apply
.. automethod:: apply_pressure
.. automethod:: apply_derivative
.. automethod:: apply_stress
"""
[docs] def __init__(self, dim=None):
self.dim = dim
if dim == 2:
self.kernel_dict = {
(3, 0): StressletKernel(dim=2, icomp=0, jcomp=0, kcomp=0),
(2, 1): StressletKernel(dim=2, icomp=0, jcomp=0, kcomp=1),
(1, 2): StressletKernel(dim=2, icomp=0, jcomp=1, kcomp=1),
(0, 3): StressletKernel(dim=2, icomp=1, jcomp=1, kcomp=1)
}
elif dim == 3:
self.kernel_dict = {
(3, 0, 0): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=0),
(2, 1, 0): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=1),
(2, 0, 1): StressletKernel(dim=3, icomp=0, jcomp=0, kcomp=2),
(1, 2, 0): StressletKernel(dim=3, icomp=0, jcomp=1, kcomp=1),
(1, 1, 1): StressletKernel(dim=3, icomp=0, jcomp=1, kcomp=2),
(1, 0, 2): StressletKernel(dim=3, icomp=0, jcomp=2, kcomp=2),
(0, 3, 0): StressletKernel(dim=3, icomp=1, jcomp=1, kcomp=1),
(0, 2, 1): StressletKernel(dim=3, icomp=1, jcomp=1, kcomp=2),
(0, 1, 2): StressletKernel(dim=3, icomp=1, jcomp=2, kcomp=2),
(0, 0, 3): StressletKernel(dim=3, icomp=2, jcomp=2, kcomp=2)
}
else:
raise ValueError("unsupported dimension given to StressletWrapper")
[docs] def apply(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit):
"""Symbolic expressions for integrating Stresslet kernel.
Returns an object array of symbolic expressions for the vector
resulting from integrating the dyadic Stresslet kernel with
variable *density_vec_sym* and source direction vectors *dir_vec_sym*.
:arg density_vec_sym: a symbolic vector variable for the density vector.
:arg dir_vec_sym: a symbolic vector variable for the direction vector.
:arg mu_sym: a symbolic variable for the viscosity.
:arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on
to :class:`~pytential.symbolic.primitives.IntG`.
"""
import itertools
sym_expr = np.empty((self.dim,), dtype=object)
for comp in range(self.dim):
# Start variable count for kernel with 1 for the requested result
# component
base_count = np.zeros(self.dim, dtype=np.int32)
base_count[comp] += 1
for i, j in itertools.product(range(self.dim), range(self.dim)):
var_ctr = base_count.copy()
var_ctr[i] += 1
var_ctr[j] += 1
ctr_key = tuple(var_ctr)
if i + j < 1:
sym_expr[comp] = sym.int_g_vec(
self.kernel_dict[ctr_key],
dir_vec_sym[i] * density_vec_sym[j],
qbx_forced_limit=qbx_forced_limit, mu=mu_sym)
else:
sym_expr[comp] = sym_expr[comp] + sym.int_g_vec(
self.kernel_dict[ctr_key],
dir_vec_sym[i] * density_vec_sym[j],
qbx_forced_limit=qbx_forced_limit,
mu=mu_sym)
return sym_expr
[docs] def apply_pressure(self, density_vec_sym, dir_vec_sym, mu_sym, qbx_forced_limit):
"""Symbolic expression for pressure field associated with the Stresslet."""
import itertools
from pytential.symbolic.mappers import DerivativeTaker
kernel = LaplaceKernel(dim=self.dim)
factor = (2. * mu_sym)
for i, j in itertools.product(range(self.dim), range(self.dim)):
if i + j < 1:
sym_expr = factor * DerivativeTaker(i).map_int_g(
DerivativeTaker(j).map_int_g(
sym.int_g_vec(kernel,
density_vec_sym[i] * dir_vec_sym[j],
qbx_forced_limit=qbx_forced_limit)))
else:
sym_expr = sym_expr + (
factor * DerivativeTaker(i).map_int_g(
DerivativeTaker(j).map_int_g(
sym.int_g_vec(kernel,
density_vec_sym[i] * dir_vec_sym[j],
qbx_forced_limit=qbx_forced_limit))))
return sym_expr
[docs] def apply_derivative(self, deriv_dir, density_vec_sym, dir_vec_sym,
mu_sym, qbx_forced_limit):
"""Symbolic derivative of velocity from stresslet.
Returns an object array of symbolic expressions for the vector
resulting from integrating the *deriv_dir* target derivative of the
dyadic Stresslet kernel with variable *density_vec_sym* and source
direction vectors *dir_vec_sym*.
:arg deriv_dir: integer denoting the axis direction for the derivative.
:arg density_vec_sym: a symbolic vector variable for the density vector.
:arg dir_vec_sym: a symbolic vector variable for the normal direction.
:arg mu_sym: a symbolic variable for the viscosity.
:arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on
to :class:`~pytential.symbolic.primitives.IntG`.
"""
import itertools
from pytential.symbolic.mappers import DerivativeTaker
sym_expr = np.empty((self.dim,), dtype=object)
for comp in range(self.dim):
# Start variable count for kernel with 1 for the requested result
# component
base_count = np.zeros(self.dim, dtype=np.int32)
base_count[comp] += 1
for i, j in itertools.product(range(self.dim), range(self.dim)):
var_ctr = base_count.copy()
var_ctr[i] += 1
var_ctr[j] += 1
ctr_key = tuple(var_ctr)
if i + j < 1:
sym_expr[comp] = DerivativeTaker(deriv_dir).map_int_g(
sym.int_g_vec(self.kernel_dict[ctr_key],
dir_vec_sym[i] * density_vec_sym[j],
qbx_forced_limit=qbx_forced_limit,
mu=mu_sym))
else:
sym_expr[comp] = sym_expr[comp] + DerivativeTaker(
deriv_dir).map_int_g(
sym.int_g_vec(self.kernel_dict[ctr_key],
dir_vec_sym[i] * density_vec_sym[j],
qbx_forced_limit=qbx_forced_limit,
mu=mu_sym))
return sym_expr
[docs] def apply_stress(self, density_vec_sym, normal_vec_sym, dir_vec_sym,
mu_sym, qbx_forced_limit):
r"""Symbolic expression for viscous stress applied to a direction.
Returns a vector of symbolic expressions for the force resulting
from the viscous stress
.. math::
-p \delta_{ij} + \mu (\nabla_i u_j + \nabla_j u_i)
applied in the direction of *dir_vec_sym*.
:arg density_vec_sym: a symbolic vector variable for the density vector.
:arg normal_vec_sym: a symbolic vector variable for the normal vectors
(outward facing normals at source locations).
:arg dir_vec_sym: a symbolic vector for the application direction.
:arg mu_sym: a symbolic variable for the viscosity.
:arg qbx_forced_limit: the *qbx_forced_limit* argument to be passed on
to :class:`~pytential.symbolic.primitives.IntG`.
"""
sym_expr = np.empty((self.dim,), dtype=object)
# Build velocity derivative matrix
sym_grad_matrix = np.empty((self.dim, self.dim), dtype=object)
for i in range(self.dim):
sym_grad_matrix[:, i] = self.apply_derivative(i, density_vec_sym,
normal_vec_sym, mu_sym, qbx_forced_limit)
for comp in range(self.dim):
# First, add the pressure term:
sym_expr[comp] = - dir_vec_sym[comp] * self.apply_pressure(
density_vec_sym, normal_vec_sym,
mu_sym, qbx_forced_limit)
# Now add the velocity derivative components
for j in range(self.dim):
sym_expr[comp] = sym_expr[comp] + (
dir_vec_sym[j] * mu_sym * (
sym_grad_matrix[comp][j]
+ sym_grad_matrix[j][comp])
)
return sym_expr
# }}}
# {{{ base Stokes operator
[docs]class StokesOperator:
"""
.. attribute:: ambient_dim
.. attribute:: side
.. automethod:: __init__
.. automethod:: get_density_var
.. automethod:: prepare_rhs
.. automethod:: operator
.. automethod:: velocity
.. automethod:: pressure
"""
[docs] def __init__(self, ambient_dim, side):
"""
:arg ambient_dim: dimension of the ambient space.
:arg side: :math:`+1` for exterior or :math:`-1` for interior.
"""
if side not in [+1, -1]:
raise ValueError(f"invalid evaluation side: {side}")
self.ambient_dim = ambient_dim
self.side = side
self.stresslet = StressletWrapper(dim=self.ambient_dim)
self.stokeslet = StokesletWrapper(dim=self.ambient_dim)
@property
def dim(self):
return self.ambient_dim - 1
[docs] def get_density_var(self, name="sigma"):
"""
:returns: a symbolic vector corresponding to the density.
"""
return sym.make_sym_vector(name, self.ambient_dim)
[docs] def prepare_rhs(self, b, *, mu):
"""
:returns: a (potentially) modified right-hand side *b* that matches
requirements of the representation.
"""
return b
[docs] def operator(self, sigma):
"""
:returns: the integral operator that should be solved to obtain the
density *sigma*.
"""
raise NotImplementedError
[docs] def velocity(self, sigma, *, normal, mu, qbx_forced_limit=None):
"""
:returns: a representation of the velocity field in the Stokes flow.
"""
raise NotImplementedError
[docs] def pressure(self, sigma, *, normal, mu, qbx_forced_limit=None):
"""
:returns: a representation of the pressure in the Stokes flow.
"""
raise NotImplementedError
# }}}
# {{{ exterior Stokes flow
[docs]class HsiaoKressExteriorStokesOperator(StokesOperator):
"""Representation for 2D Stokes Flow based on [HsiaoKress1985]_.
Inherits from :class:`StokesOperator`.
.. [HsiaoKress1985] G. C. Hsiao and R. Kress, *On an Integral Equation for
the Two-Dimensional Exterior Stokes Problem*,
Applied Numerical Mathematics, Vol. 1, 1985,
`DOI <https://doi.org/10.1016/0168-9274(85)90029-7>`__.
.. automethod:: __init__
"""
[docs] def __init__(self, *, omega, alpha=None, eta=None):
r"""
:arg omega: farfield behaviour of the velocity field, as defined
by :math:`A` in [HsiaoKress1985]_ Equation 2.3.
:arg alpha: real parameter :math:`\alpha > 0`.
:arg eta: real parameter :math:`\eta > 0`. Choosing this parameter well
can have a non-trivial effect on the conditioning.
"""
super().__init__(ambient_dim=2, side=+1)
# NOTE: in [hsiao-kress], there is an analysis on a circle, which
# recommends values in
# 1/2 <= alpha <= 2 and max(1/alpha, 1) <= eta <= min(2, 2/alpha)
# so we choose alpha = eta = 1, which seems to be in line with some
# of the presented numerical results too.
if alpha is None:
alpha = 1.0
if eta is None:
eta = 1.0
self.omega = omega
self.alpha = alpha
self.eta = eta
def _farfield(self, mu, qbx_forced_limit):
length = sym.integral(self.ambient_dim, self.dim, 1)
return self.stokeslet.apply(
-self.omega / length,
mu,
qbx_forced_limit=qbx_forced_limit)
def _operator(self, sigma, normal, mu, qbx_forced_limit):
slp_qbx_forced_limit = qbx_forced_limit
if slp_qbx_forced_limit == "avg":
slp_qbx_forced_limit = +1
# NOTE: we set a dofdesc here to force the evaluation of this integral
# on the source instead of the target when using automatic tagging
# see :meth:`pytential.symbolic.mappers.LocationTagger._default_dofdesc`
dd = sym.DOFDescriptor(None, discr_stage=sym.QBX_SOURCE_STAGE1)
int_sigma = sym.integral(self.ambient_dim, self.dim, sigma, dofdesc=dd)
meanless_sigma = sym.cse(sigma - sym.mean(self.ambient_dim, self.dim, sigma))
op_k = self.stresslet.apply(sigma, normal, mu,
qbx_forced_limit=qbx_forced_limit)
op_s = (
self.alpha / (2.0 * np.pi) * int_sigma
- self.stokeslet.apply(meanless_sigma, mu,
qbx_forced_limit=slp_qbx_forced_limit)
)
return op_k + self.eta * op_s
def prepare_rhs(self, b, *, mu):
return b + self._farfield(mu, qbx_forced_limit=+1)
def operator(self, sigma, *, normal, mu):
# NOTE: H. K. 1985 Equation 2.18
return -0.5 * self.side * sigma - self._operator(sigma, normal, mu, "avg")
def velocity(self, sigma, *, normal, mu, qbx_forced_limit=2):
# NOTE: H. K. 1985 Equation 2.16
return (
-self._farfield(mu, qbx_forced_limit)
- self._operator(sigma, normal, mu, qbx_forced_limit)
)
def pressure(self, sigma, *, normal, mu, qbx_forced_limit=2):
# FIXME: H. K. 1985 Equation 2.17
raise NotImplementedError
[docs]class HebekerExteriorStokesOperator(StokesOperator):
"""Representation for 3D Stokes Flow based on [Hebeker1986]_.
Inherits from :class:`StokesOperator`.
.. [Hebeker1986] F. C. Hebeker, *Efficient Boundary Element Methods for
Three-Dimensional Exterior Viscous Flow*, Numerical Methods for
Partial Differential Equations, Vol. 2, 1986,
`DOI <https://doi.org/10.1002/num.1690020404>`__.
.. automethod:: __init__
"""
[docs] def __init__(self, *, eta=None):
r"""
:arg eta: a parameter :math:`\eta > 0`. Choosing this parameter well
can have a non-trivial effect on the conditioning of the operator.
"""
super().__init__(ambient_dim=3, side=+1)
# NOTE: eta is chosen here based on H. 1986 Figure 1, which is
# based on solving on the unit sphere
if eta is None:
eta = 0.75
self.eta = eta
def _operator(self, sigma, normal, mu, qbx_forced_limit):
slp_qbx_forced_limit = qbx_forced_limit
if slp_qbx_forced_limit == "avg":
slp_qbx_forced_limit = self.side
op_w = self.stresslet.apply(sigma, normal, mu,
qbx_forced_limit=qbx_forced_limit)
op_v = self.stokeslet.apply(sigma, mu,
qbx_forced_limit=slp_qbx_forced_limit)
return op_w + self.eta * op_v
def operator(self, sigma, *, normal, mu):
# NOTE: H. 1986 Equation 17
return -0.5 * self.side * sigma - self._operator(sigma, normal, mu, "avg")
def velocity(self, sigma, *, normal, mu, qbx_forced_limit=2):
# NOTE: H. 1986 Equation 16
return -self._operator(sigma, normal, mu, qbx_forced_limit)
def pressure(self, sigma, *, normal, mu, qbx_forced_limit=2):
# FIXME: not given in H. 1986, but should be easy to derive using the
# equivalent single-/double-layer pressure kernels
raise NotImplementedError
# }}}