yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #22884
Re: [Question #690127]: Missing particles in VTK
Question #690127 on Yade changed:
https://answers.launchpad.net/yade/+question/690127
Przemek posted a new comment:
import math
import numpy as np
from scipy.sparse.linalg import spsolve
from numpy.linalg import norm
from numpy.linalg import inv
from yade import pack, export
# Properties and parameters
kappa = 10. # heat conductivity
beta = 1.
T_inf = 20. # ambient temperature
q_v = 0.5e1 # heat flux
convection1 = 0. # x-position of convection
convection2 = 0. # y-position of convection
tPoint = [2, 2, 0]
avTemp = 0
# Plate dimension
boxMax = (5.0, 5.0, 0.5)
# Spheres
color = [0.21, 0.22, 0.1]
sphereRadius = 0.1 # cm
r = 0.11
ndiv_x = boxMax[0] / (2 * sphereRadius) + 1
ndiv_y = boxMax[1] / (2 * sphereRadius) + 1
ndiv_z = boxMax[2] / (2 * sphereRadius) + 1
nbSpheres = (int(ndiv_x), int(ndiv_y), int(ndiv_z))
num = 0
for i in range(nbSpheres[0]):
for j in range(nbSpheres[1]):
for jj in range(nbSpheres[2]):
x = (2.0 * sphereRadius * i)
y = (2.0 * sphereRadius * j)
z = (2.0 * sphereRadius * jj)
num += 1
plate = utils.sphere([x, y, z], r, color=color, fixed=fixed)
O.bodies.append(plate)
conv1 = box(center=(convection1, boxMax[1] / 2, boxMax[2] / 2), extents=(0., boxMax[1] / 2 + sphereRadius, boxMax[2] / 2 + sphereRadius),
fixed=True, wire=True, mask=3),
conv2 = box(center=(boxMax[0] / 2, convection2, boxMax[2] / 2), extents=(boxMax[0] / 2 + sphereRadius, 0., boxMax[2] / 2 + sphereRadius),
fixed=True, wire=True, mask=3),
heat1 = facet(vertices=[(boxMax[0] + sphereRadius, -sphereRadius, -sphereRadius),
(boxMax[0] + sphereRadius, boxMax[0] + sphereRadius, -sphereRadius),
(boxMax[0] + sphereRadius, -sphereRadius, boxMax[2] + sphereRadius)], fixed=True, wire=True, mask=55)
heat2 = facet(vertices=[(boxMax[0] + sphereRadius, boxMax[0] + sphereRadius, boxMax[2] + sphereRadius),
(boxMax[0] + sphereRadius, boxMax[0] + sphereRadius, -sphereRadius),
(boxMax[0] + sphereRadius, -sphereRadius, boxMax[2] + sphereRadius)], fixed=True, wire=True, mask=55)
O.bodies.append(heat1)
O.bodies.append(heat2)
O.bodies.append(conv1)
O.bodies.append(conv2)
nBox = 2 # number of boxes
###
def convComp():
betaDiffx = 0
betaDiffy = 0
for b in O.bodies:
if isinstance(b.shape, Sphere):
if b.state.pos[0] == convection1:
betaDiffx += 1
elif b.state.pos[1] == convection2:
betaDiffy += 1
else:
pass
beta_x = beta / betaDiffx
beta_y = beta / betaDiffy
return beta_x, beta_y
###
beta_x = convComp()[0]
beta_y = convComp()[1]
# Solver
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb(), Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
NewtonIntegrator(gravity=(0, 0, -9.81), damping=0.1),
PyRunner(command='intCounter()', realPeriod=1),
PyRunner(command='tempCalc()', realPeriod=1),
PyRunner(command='checkTempPoint(tPoint, sphereRadius)', realPeriod=1)
]
# Make a step
O.step()
###
def intCounter(): # counter of body (Sphere and Facet) real interactions
for b in O.bodies:
if isinstance(b.shape, Sphere) == True:
countI = 0
bInCon = []
for i in O.interactions:
if i.isReal:
if (i.id1 == b.id):
if isinstance(O.bodies[i.id2].shape, Facet) == True:
countI += 1
bInCon.append(i.id2)
else:
pass
elif (i.id2 == b.id):
if isinstance(O.bodies[i.id1].shape, Facet) == True:
countI += 1
bInCon.append(i.id1)
else:
pass
else:
pass
b.realInt(countI)
b.bodiesInContact(bInCon)
elif isinstance(b.shape, Facet) == True:
countI = 0
bInCon = []
for i in O.interactions:
if i.isReal:
if (i.id1 == b.id):
countI += 1
bInCon.append(i.id2)
elif (i.id2 == b.id):
countI += 1
bInCon.append(i.id1)
else:
pass
b.realInt(countI)
b.bodiesInContact(bInCon)
else:
pass
###
###
def volumefraction(fID, sID):
ns = O.bodies[fID].countI
nsd = 0
# Ai = pi * math.pow(r * math.sin(math.acos((r - O.interactions[sID, fID].geom.penetrationDepth) / r)), 2) # Area of penetraded particle
for b in O.bodies[fID].bInCon:
if O.bodies[b].countI > 1:
nsd += 1
else:
pass
nsn = ns - nsd
if O.bodies[sID].countI > 1:
vf = [O.interactions[sID, fID].geom.penetrationDepth]
for fCon in O.bodies[sID].bInCon:
if fCon != fID:
vf.append(O.interactions[sID, fCon].geom.penetrationDepth)
else:
pass
psi = vf[0] / sum(vf) / ns
else:
psi_prime = []
for sCon in O.bodies[fID].bInCon:
if sCon != sID:
if O.bodies[sCon].countI > 1:
vf = [O.interactions[sCon, fID].geom.penetrationDepth]
for fCon in O.bodies[sCon].bInCon:
if fCon != fID:
vf.append(O.interactions[sCon, fCon].geom.penetrationDepth)
else:
pass
psi_prime.append(vf[0] / sum(vf))
else:
pass
else:
pass
psi = (1 - (sum(psi_prime) / ns)) / float(nsn)
return psi
###
###
def tempCalc():
def interNumber(): # is needed???
tmp = 0
vec = []
for i in O.interactions:
if i.isReal:
vec.append([i.id1, i.id2])
tmp += 1
else:
pass
return tmp, vec
realInt = interNumber()[0] # real number of interactions
interList = interNumber()[1] # list of bodies IDs in interactions
def interactionsDef(ii, jj):
conPoint = O.interactions[ii, jj].geom.contactPoint # coordinates of contact point
penDist = O.interactions[ii, jj].geom.penetrationDepth / 2. # penetration distance of current particle
alphaI = math.acos(1 - penDist / r) # contact angle - alpha_i
###
# ui = (conPoint - xi) / 2.
xi = O.bodies[ii].state.pos # current position of body 1
xj = O.bodies[jj].state.pos # current position of body 1
ri = xj - xi
ui = np.array([1., 0., 0.]) # basis vector x
vi = np.array([0., 1., 0.]) # basis vector y
wi = np.array([0., 0., 1.]) # basis vector z
teta_1 = np.arccos(
(ri.dot(vi.T)) / (norm(ri) * norm(vi))) # angle of contact position theta_i - direction 1
teta_2 = np.arccos(
(ri.dot(wi.T)) / (norm(ri) * norm(wi))) # angle of contact position theta_i - direction 2
cx = (ri.dot(ui.T)) / (norm(ri) * norm(ui))
cy = (ri.dot(vi.T)) / (norm(ri) * norm(vi))
cz = (ri.dot(wi.T)) / (norm(ri) * norm(wi))
c11 = cx * cx
c22 = cy * cy
c33 = cz * cz
c12 = cx * cy
c13 = cx * cz
c23 = cy * cz
cosMatrix = np.array([[c11, c12, c13, -c11, -c12, -c13],
[c12, c22, c23, -c12, -c22, -c23],
[c13, c23, c33, -c13, -c23, -c33],
[-c11, -c12, -c13, c11, c12, c13],
[-c12, -c22, -c23, c12, c22, c23],
[-c13, -c23, -c33, c13, c23, c33]])
return alphaI, cosMatrix
H = np.zeros([len(O.bodies) - nBox, len(O.bodies) - nBox]) # definition of thermal resistance matrix
F = np.zeros(len(O.bodies) - nBox) # definition of thermal force vector
def thResMatrix():
for con in O.interactions: # loop over all interactions
if con.isReal:
idBodies = [con.id1, con.id2] # Allocation vector
ii = idBodies[0] # ID of body 1
jj = idBodies[1] # ID of body 2
if isinstance(O.bodies[ii].shape, Sphere) and isinstance(O.bodies[jj].shape,
Sphere) == True: # only spheres
alphaI = interactionsDef(ii, jj)[0]
K_el = np.power(2 / math.pi / kappa * (3 / 2 + np.power(alphaI, 2) / 36.
+ np.power(alphaI, 4) / 2700. + np.power(alphaI, 6) / 79380.
+ np.power(alphaI, 8) / 252000. - np.log(alphaI)), -1)
H[ii][ii] += K_el
H[ii][jj] += -K_el
H[jj][ii] += -K_el
H[jj][jj] += K_el
for conv in range(2):
if O.bodies[idBodies[conv]].state.pos[0] == convection1:
H[idBodies[conv]][idBodies[conv]] += beta_x
F[idBodies[conv]] += beta_x * T_inf
if O.bodies[idBodies[conv]].state.pos[1] == convection2:
H[idBodies[conv]][idBodies[conv]] += beta_y
F[idBodies[conv]] += beta_y * T_inf
else:
pass
elif O.bodies[idBodies[conv]].state.pos[1] == convection2:
H[idBodies[conv]][idBodies[conv]] += beta_y
F[idBodies[conv]] += beta_y * T_inf
if O.bodies[idBodies[conv]].state.pos[0] == convection1:
H[idBodies[conv]][idBodies[conv]] += beta_x
F[idBodies[conv]] += beta_x * T_inf
else:
pass
else:
pass
elif (isinstance(O.bodies[ii].shape, Facet) or isinstance(O.bodies[jj].shape, Facet)) == True: # facets
alphaI = interactionsDef(ii, jj)[0]
K_el = np.power(1 / math.pi / kappa * (3 / 2 + np.power(alphaI, 2) / 36.
+ np.power(alphaI, 4) / 2700. + np.power(alphaI, 6) / 79380.
+ np.power(alphaI, 8) / 252000. - np.log(alphaI)), -1)
H[ii][ii] += K_el
H[ii][jj] += -K_el
H[jj][ii] += -K_el
H[jj][jj] += K_el
for conv in range(2):
if isinstance(O.bodies[idBodies[conv]].shape, Facet) == True:
volFrac = volumefraction(idBodies[conv], idBodies[1 - conv])
F[idBodies[conv]] += q_v * volFrac
else:
pass
else:
pass
else:
pass
def boundaryCon():
for bon in O.interactions: # loop over all bodies
if bon.isReal:
idBodies = [bon.id1, bon.id2] # Allocation vector
if isinstance(O.bodies[idBodies[0]].shape, Sphere) or isinstance(O.bodies[idBodies[1]].shape,
Sphere) != True:
ii = i.id # body ID
xi = i.state.pos # current position of body
for j in range(len(i.intrs())): # loop over all bodies in interaction with current body
k = utils.getBodyIdsContacts(ii)[j] # ID of contact particle
if O.interactions[ii, k].isReal:
if isinstance(i.shape, Sphere): # only spheres
alphaI = interactionsDef(ii, xi, k)
H[ii] += 2 / math.pi / kappa * (3 / 2 + np.power(alphaI, 2) / 36.
+ np.power(alphaI, 4) / 2700. + np.power(alphaI,
6) / 79380.
+ np.power(alphaI, 8) / 252000. - np.log(alphaI))
else: # other shapes - facets, boxes
alphaI = interactionsDef(ii, xi, k)
H[k] += 1 / math.pi / kappa * (3 / 2 + np.power(alphaI, 2) / 36.
+ np.power(alphaI, 4) / 2700. + np.power(alphaI,
6) / 79380.
+ np.power(alphaI, 8) / 252000. - np.log(alphaI))
if i.mask == 45:
H[k] += + beta
else:
pass
else:
pass
thResMatrix()
#T = inv(H).dot(F) # add a module to solve the sparse matrix
T = spsolve(H, F)
minT = min(T)
maxT = max(T)
#ran = floor(max(T) - minT)
ran = maxT - minT
for tp in O.bodies:
if isinstance(tp.shape, Box) == True:
pass
else:
tp.setTemp(T[tp.id])
colMap = (T[tp.id] - minT) / ran * 255
n = 4 * colMap / 256
red = 255 * min(max(min(n - 1.5, -n + 4.5), 0), 1)
green = 255 * min(max(min(n - 0.5, -n + 3.5), 0), 1)
blue = 255 * min(max(min(n + 0.5, -n + 2.5), 0), 1)
tp.shape.color = (red, green, blue)
export.textExt('test.txt', format='x_y_z_r_attrs', attrs=['b.getTemp()'], comment='temp')
export.text2vtk('test.txt', 'test.vtk')
print(maxT - minT)
###
###
def checkTempPoint(tPoint, r):
scaleFactor = 1.5
inflRad = scaleFactor * r
avTemp = []
for b in O.bodies:
if type(b.shape) == Sphere:
bRad = norm(b.state.pos - tPoint)
if bRad <= inflRad:
avTemp.append(b.getTemp())
else:
pass
avTemp = np.mean(avTemp)
print(avTemp)
###
from yade import pack, qt
v = qt.View()
v.axes = False
--
You received this question notification because your team yade-users is
an answer contact for Yade.
