In [1]:

    

# !/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on 

@author: zhangji
"""

%pylab inline
pylab.rcParams['figure.figsize'] = (18.5, 10.5)
fontsize = 40

import os
import importlib
from time import time
import numpy as np
import scipy as sp
import pandas as pd
import re
from scanf import scanf
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d, Axes3D
from scipy.optimize import leastsq, curve_fit
from IPython.display import display, HTML
from scipy import interpolate, integrate, optimize, sparse
from codeStore import support_fun as spf
from src import slenderBodyTheory as slb

PWD = os.getcwd()
np.set_printoptions(linewidth=130, precision=5)


    

    




Populating the interactive namespace from numpy and matplotlib

In [74]:

    

importlib.reload(slb)

ph = 1
rt1 = 1
rt2 = 1
theta = 0
uz, wz = 0, 1

phi_list = np.linspace(0, 2 * np.pi, 30)
t_fra = slb.T1_fun(theta, ph, rt1, rt2)
n_fra = slb.N1_fun(theta, ph, rt1, rt2)
b_fra = slb.B1_fun(theta, ph, rt1, rt2)
ep = np.array([n_fra * ti for ti in np.cos(phi_list)]) + np.array([b_fra * ti for ti in np.sin(phi_list)])
x0 = slb.x1_fun(theta, ph, rt1, rt2)
x1 = x0 + ep * rt2
u1 = np.array((0, 0, uz)) + np.cross(np.array((0, 0, wz)), x1)
fft_u = np.fft.fft(u1, axis=0)
usin = -fft_u.imag
ucos = fft_u.real
# Eq 5.4 in Koens2018
t1 = 2 * (np.eye(3) + np.outer(t_fra, t_fra)) - np.outer(n_fra, n_fra) + np.outer(b_fra, b_fra)
t2 = -(np.outer(n_fra, b_fra) + np.outer(b_fra, n_fra))
t3 = 2 * (np.eye(3) + np.outer(t_fra, t_fra)) + np.outer(n_fra, n_fra) - np.outer(b_fra, b_fra)
t4 = -(np.outer(n_fra, b_fra) + np.outer(b_fra, n_fra))
tA = np.vstack((np.hstack((t1, t2)), np.hstack((t3, t4))))
tb = np.hstack((8 * ucos[1], 8 * usin[1]))


    

In [9]:

    

importlib.reload(slb)
x1_fun = slb.x1_fun
T1_fun = slb.T1_fun

def KRJ_stokeslets_mij2(u_theta, f_theta, fidx, ph, rt1, rt2,
                        u_node_fun=x1_fun, f_node_fun=x1_fun, T_fun=T1_fun):
    # inner_mj = S(:, j), along u
    S = np.sqrt(4 * np.pi ** 2 * rt1 ** 2 + ph ** 2)
    intFct = S / (2 * np.pi)
    su = u_theta * intFct
    sf = f_theta * intFct
    t = T_fun(u_theta, ph, rt1, rt2)
    ds = np.abs(sf - su)
    if u_node_fun is f_node_fun:
        ds[fidx] = np.inf
    t_m = np.vstack([(np.eye(3) + np.outer(ti, ti)) / dsi for ti, dsi in zip(t.reshape(-1, 3), ds)])
    return t_m

ph = 1
rt1 = 1
rt2 = 1
u_theta = np.linspace(0, 2 * np.pi, 3)
# f_theta = 0.3 * np.pi
fidx=1

tm = KRJ_stokeslets_mij2(u_theta, u_theta[fidx], fidx, ph, rt1, rt2, 
                         u_node_fun=x1_fun, f_node_fun=x1_fun, T_fun=T1_fun)


    

    




[[ 3.14353e-01  0.00000e+00  0.00000e+00]
 [ 0.00000e+00  5.02659e+02  7.99506e+01]
 [ 0.00000e+00  7.99506e+01  1.30389e+01]
 [ 0.00000e+00  0.00000e+00 -0.00000e+00]
 [ 0.00000e+00  0.00000e+00 -0.00000e+00]
 [-0.00000e+00 -0.00000e+00  0.00000e+00]
 [ 3.14353e-01  1.23039e-13  1.95822e-14]
 [ 1.23039e-13  5.02659e+02  7.99506e+01]
 [ 1.95822e-14  7.99506e+01  1.30389e+01]]

In [60]:

    

Fn1Mat_fun = lambda theta, ph, rt1, rt2: \
    np.identity(3) - np.tensordot(T1_fun(theta, ph, rt1, rt2), T1_fun(theta, ph, rt1, rt2), axes=0)
Fn2Mat_fun = lambda theta, ph, rt1, rt2: \
    np.identity(3) - np.tensordot(T2_fun(theta, ph, rt1, rt2), T2_fun(theta, ph, rt1, rt2), axes=0)


    

    
Out[60]:





array([ 0.25016, 14.63443,  0.     ])

In [38]:

    

T_fun=T1_fun
u_theta = np.array((0.1, 0.2))
t = T_fun(u_theta, ph, rt1, rt2)
print(t)
print()
temp1 = np.vstack([np.eye(3) + np.outer(ti, ti) for ti in t.reshape(-1, 3)])
print(temp1)


    

    




[[-0.09859  0.98264  0.15718]
 [-0.1962   0.96788  0.15718]]

[[ 1.00972 -0.09688 -0.0155 ]
 [-0.09688  1.96557  0.15445]
 [-0.0155   0.15445  1.0247 ]
 [ 1.03849 -0.1899  -0.03084]
 [-0.1899   1.9368   0.15213]
 [-0.03084  0.15213  1.0247 ]]

In [ ]:

    

[[-0.09859  0.98264  0.15718]]

[[ 0.00972 -0.09688 -0.0155 ]
 [-0.09688  0.96557  0.15445]
 [-0.0155   0.15445  0.0247 ]]


    

In [ ]:

    

[[-0.1962   0.96788  0.15718]]

[[ 0.03849 -0.1899  -0.03084]
 [-0.1899   0.9368   0.15213]
 [-0.03084  0.15213  0.0247 ]]