/usr/lib/python2.7/dist-packages/PySPH-1.0a4.dev0-py2.7-linux-x86_64.egg/pysph/sph/boundary_equations.py is in python-pysph 0~20160514.git91867dc-4build1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 | """
SPH Boundary Equations
######################
"""
from pysph.sph.equation import Equation
def wendland_quintic(rij=1.0, h=1.0):
q = rij/h
q1 = 2.0 - q
val = 0.0
if q < 2.0:
val = (1 + 2.5*q + 2*q*q)*q1*q1*q1*q1*q1
return val
class MonaghanBoundaryForce(Equation):
def __init__(self, dest, sources, deltap):
self.deltap = deltap
super(MonaghanBoundaryForce,self).__init__(dest,sources)
def loop(self, d_idx, s_idx, s_m, s_rho, d_m, d_cs, s_cs, d_h,
s_tx, s_ty, s_tz, s_nx, s_ny, s_nz,
d_au, d_av, d_aw, XIJ):
norm = declare('matrix((3,))')
tang = declare('matrix((3,))')
ma = d_m[d_idx]
mb = s_m[s_idx]
# particle sound speed
cs = d_cs[d_idx]
# boundary normals
norm[0] = s_nx[s_idx]
norm[1] = s_ny[s_idx]
norm[2] = s_nz[s_idx]
# boundary tangents
tang[0] = s_tx[s_idx]
tang[1] = s_ty[s_idx]
tang[2] = s_tz[s_idx]
# x and y projections
x = XIJ[0]*tang[0] + XIJ[1]*tang[1] + XIJ[2]*tang[2]
y = XIJ[0]*norm[0] + XIJ[1]*norm[1] + XIJ[2]*norm[2]
# compute the force
force = 0.0
q = y/d_h[d_idx]
xabs = fabs(x)
if 0 <= xabs <= self.deltap:
beta = 0.02 * cs * cs/y
tforce = 1.0 - xabs/self.deltap
if 0 < q <= 2.0/3.0:
nforce = 2.0/3.0
elif 2.0/3.0 < q <= 1.0:
nforce = 2*q*(1.0 - 0.75*q)
elif 1.0 < q <= 2.0:
nforce = 0.5 * (2-q)*(2-q)
else:
nforce = 0.0
force = (mb/(ma+mb)) * nforce * tforce * beta
else:
force = 0.0
# boundary force accelerations
d_au[d_idx] += force * norm[0]
d_av[d_idx] += force * norm[1]
d_aw[d_idx] += force * norm[2]
class MonaghanKajtarBoundaryForce(Equation):
def __init__(self, dest, sources, K=None, beta=None, h=None):
self.K = K
self.beta = beta
self.h = h
if None in [K, beta, h]:
raise ValueError("Invalid parameter values")
super(MonaghanKajtarBoundaryForce,self).__init__(dest,sources)
def _get_helpers_(self):
return [wendland_quintic]
def loop(self, d_idx, s_idx, d_m, s_m, d_au, d_av, d_aw, RIJ, R2IJ, XIJ):
ma = d_m[d_idx]
mb = s_m[s_idx]
w = wendland_quintic(RIJ, self.h)
force = self.K/self.beta * w/R2IJ * 2*mb/(ma + mb)
d_au[d_idx] += force * XIJ[0]
d_av[d_idx] += force * XIJ[1]
d_aw[d_idx] += force * XIJ[2]
|