added hydrodynamic interface code

stoked/__init__.py

@@ -14,7 +14,7 @@
 from .electrodynamics import point_dipole_electrodynamics
 from .van_der_waals import van_der_waals_sphere, van_der_waals_sphere_interface
 from .gravity import gravity, gravity_sphere, gravity_ellipsoid
-from .hydrodynamics import grand_mobility_matrix
+from .hydrodynamics import grand_mobility_matrix, interface
 from .constraints import constrain_position, constrain_rotation
 from .forces import pairwise_central_force, pairwise_force
 from .common import lennard_jones

stoked/hydrodynamics.py

@@ -1,5 +1,14 @@
 import numpy as np
 
+class interface:
+    """A no-slip interface"""
+    def __init__(self, z=0):
+        """
+        Arguments:
+            z   z-position of the interface
+        """
+        self.z = z
+
 def levi_civita():
     """return the levi-civita symbol"""
 
@@ -8,6 +17,32 @@ def levi_civita():
     eijk[0, 2, 1] = eijk[2, 1, 0] = eijk[1, 0, 2] = -1
     return eijk
 
+def particle_wall_self_mobility(position, interface, viscosity, radius):
+    """
+    Construct the particle wall self-mobility matrix for a single particle
+
+    Arguments:
+        position[3]       position of particle
+        interface         interface object
+        viscosity         dynamic viscosity µ of surrounding fluid
+        radius            particle radius
+    """
+    M = np.zeros([2, 2, 3, 3], dtype=float)
+    h = position[2] - interface.z
+
+    gamma_T = 6*np.pi*viscosity*radius
+    gamma_R = 6*np.pi*viscosity*radius**3
+
+    a = 1/(16*gamma_T)*(9/h - 2/h**3 + 1/h**5)
+    b = 1/(8*gamma_T)*(9/h - 4/h**3 + 1/h**5)
+    M[0,0] = np.diag([a,a,b])
+
+    a = 15/(64*gamma_R)*(1/h**3)
+    b = 3/(32*gamma_R)*(1/h**3)
+    M[1,1] = np.diag([a,a,b])
+
+    return M
+
 def grand_mobility_matrix(position, drag_T, drag_R, viscosity):
     """
     Construct the grand mobility matrix for a given cluster

stoked/solver.py

@@ -63,7 +63,7 @@ class stokesian_dynamics:
     """
     def __init__(self, *, temperature, dt, position, drag, orientation=None, 
                  force=None, torque=None, interactions=None, constraint=None,
-                 hydrodynamic_coupling=True):
+                 interface=None, hydrodynamic_coupling=True):
         """
         Arguments:
             temperature        system temperature
@@ -75,6 +75,7 @@ def __init__(self, *, temperature, dt, position, drag, orientation=None,
             torque(t, r, q)    external torque function given time t, position r[N,D], orientation q[N] and returns torque T[N,D] (can be a list of functions)
             interactions       particle interactions (can be a list)
             constraint         particle constraints (can be a list)
+            interface          no-slip boundary interface (default: no interface)
             hydrodynamic_coupling    if True, include hydrodynamic coupling interactions (default: True)
         """
         self.temperature = temperature
@@ -84,6 +85,7 @@ def __init__(self, *, temperature, dt, position, drag, orientation=None,
         self.drag = drag
         self.Nparticles = len(self.position)
         self.ndim = self.position.shape[1]
+        self.interface = interface
         self.hydrodynamic_coupling = hydrodynamic_coupling
 
         if orientation is not None: