# encoding: utf-8
# 2008 © Václav Šmilauer <eudoxos@arcig.cz>
#
# I doubt there functions will be useful for anyone besides me.
#
"""Miscillaneous functions that are not believed to be generally usable,
therefore kept in my "private" module here.
They comprise notably oofem export and various CPM-related functions.
"""
from yade.wrapper import *
from math import *
from yade._eudoxos import * ## c++ implementations
[docs]class IntrSmooth3d():
r"""Return spatially weigted gaussian average of arbitrary quantity defined on interactions.
At construction time, all real interactions are put inside spatial grid, permitting fast search for
points in neighbourhood defined by distance.
Parameters for the distribution are standard deviation :math:`\sigma` and relative cutoff distance
*relThreshold* (3 by default) which will discard points farther than *relThreshold* :math:`\times \sigma`.
Given central point :math:`p_0`, points are weighted by gaussian function
.. math::
\rho(p_0,p)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(\frac{-||p_0-p||^2}{2\sigma^2}\right)
To get the averaged value, simply call the instance, passing central point and callable object which received interaction object and returns the desired quantity:
>>> O.reset()
>>> from yade import utils
>>> O.bodies.append([utils.sphere((0,0,0),1),utils.sphere((0,0,1.9),1)])
[0, 1]
>>> O.engines=[InteractionLoop([Ig2_Sphere_Sphere_Dem3DofGeom(),],[Ip2_FrictMat_FrictMat_FrictPhys()],[])]
>>> utils.createInteraction(0,1) #doctest: +ELLIPSIS
<Interaction instance at 0x...>
>> is3d=IntrSmooth3d(0.003)
>> is3d((0,0,0),lambda i: i.phys.normalForce)
Vector3(0,0,0)
"""
def __init__(self,stDev):
self.locator=InteractionLocator()
self.stDev=stDev
self.relThreshold=3.
self.sqrt2pi=sqrt(2*pi)
import yade.config
if not 'vtk' in yade.config.features: raise RuntimeError("IntrSmooth3d is only function with VTK-enabled builds.")
def _ptpt2weight(self,pt0,pt1):
distSq=(pt0-pt1).SquaredLength()
return (1./(self.stDev*self.sqrt2pi))*exp(-distSq/(2*self.stDev*self.stDev))
[docs] def bounds(self): return self.locator.bounds()
[docs] def count(self): return self.locator.count()
def __call__(self,pt,extr):
intrs=self.locator.intrsAroundPt(pt,self.stDev*self.relThreshold)
if len(intrs)==0: return float('nan')
weights,val=0.,0.
for i in intrs:
weight=self._ptpt2weight(pt,i.geom.contactPoint)
val+=weight*extr(i)
weights+=weight
#print i,weight,extr(i)
return val/weights
[docs]def estimateStress(strain,cutoff=0.):
"""Use summed stored energy in contacts to compute macroscopic stress over the same volume, provided known strain."""
# E=(1/2)σεAl # global stored energy
# σ=EE/(.5εAl)=EE/(.5εV)
from yade import utils
dim=utils.aabbDim(cutoff,centers=False)
return utils.elasticEnergy(utils.aabbExtrema(cutoff))/(.5*strain*dim[0]*dim[1]*dim[2])
[docs]def estimatePoissonYoung(principalAxis,stress=0,plot=False,cutoff=0.):
"""Estimate Poisson's ration given the "principal" axis of straining.
For every base direction, homogenized strain is computed
(slope in linear regression on discrete function particle coordinate →
→ particle displacement in the same direction as returned by
utils.coordsAndDisplacements) and, (if axis '0' is the strained
axis) the poisson's ratio is given as -½(ε₁+ε₂)/ε₀.
Young's modulus is computed as σ/ε₀; if stress σ is not given (default 0),
the result is 0.
cutoff, if > 0., will take only smaller part (centered) or the specimen into account
"""
dd=[] # storage for linear regression parameters
import pylab,numpy
try:
import stats
except ImportError:
raise ImportError("Unable to import stats; install the python-stats package.")
from yade import utils
if cutoff>0: cut=utils.fractionalBox(fraction=1-cutoff)
for axis in [0,1,2]:
if cutoff>0:
w,dw=utils.coordsAndDisplacements(axis,Aabb=cut)
else:
w,dw=utils.coordsAndDisplacements(axis)
l,ll=stats.linregress(w,dw)[0:2] # use only tangent and section
dd.append((l,ll,min(w),max(w)))
if plot: pylab.plot(w,dw,'.',label=r'$\Delta %s(%s)$'%('xyz'[axis],'xyz'[axis]))
if plot:
for axis in [0,1,2]:
dist=dd[axis][-1]-dd[axis][-2]
c=numpy.linspace(dd[axis][-2]-.2*dist,dd[axis][-1]+.2*dist)
d=[dd[axis][0]*cc+dd[axis][1] for cc in c]
pylab.plot(c,d,label=r'$\widehat{\Delta %s}(%s)$'%('xyz'[axis],'xyz'[axis]))
pylab.legend(loc='upper left')
pylab.xlabel(r'$x,\;y,\;z$')
pylab.ylabel(r'$\Delta x,\;\Delta y,\; \Delta z$')
pylab.show()
otherAxes=(principalAxis+1)%3,(principalAxis+2)%3
avgTransHomogenizedStrain=.5*(dd[otherAxes[0]][0]+dd[otherAxes[1]][0])
principalHomogenizedStrain=dd[principalAxis][0]
return -avgTransHomogenizedStrain/principalHomogenizedStrain,stress/principalHomogenizedStrain
[docs]def oofemTextExport(fName):
"""Export simulation data in text format
The format is line-oriented as follows::
E G # elastic material parameters
epsCrackOnset crackOpening xiShear transStrainCoeff # tensile parameters; epsFr=crackOpening/len
cohesionT tanPhi # shear parameters
number_of_spheres number_of_links
id x y z r boundary # spheres; boundary: -1 negative, 0 none, 1 positive
…
id1 id2 cp_x cp_y cp_z A # interactions; cp = contact point; A = cross-section
"""
material,bodies,interactions=[],[],[]
f=open(fName,'w') # fail early on I/O problem
ph=O.interactions.nth(0).phys # some params are the same everywhere
material.append("%g %g"%(ph.E,ph.G))
material.append("%g %g %g %g"%(ph.epsCrackOnset,ph.epsFracture,1e50,0.0))
material.append("%g %g"%(ph.undamagedCohesion,ph.tanFrictionAngle))
# need strainer for getting bodies in positive/negative boundary
strainers=[e for e in O.engines if e.name=='UniaxialStrainer']
if len(strainers)>0: strainer=strainers[0]
else: strainer=None
for b in O.bodies:
if not b.shape or b.shape.name!='Sphere': continue
if strainer and b.id in strainer.negIds: boundary=-1
elif strainer and b.id in strainer.posIds: boundary=1
else: boundary=0
bodies.append("%d %g %g %g %g %d"%(b.id,b.state.pos[0],b.state.pos[1],b.state.pos[2],b.shape.radius,boundary))
for i in O.interactions:
cp=i.geom.contactPoint
interactions.append("%d %d %g %g %g %g"%(i.id1,i.id2,cp[0],cp[1],cp[2],i.phys.crossSection))
f.write('\n'.join(material+["%d %d"%(len(bodies),len(interactions))]+bodies+interactions))
f.close()
[docs]def oofemPrescribedDisplacementsExport(fileName):
f=open(fileName,'w')
f.write(fileName+'.out\n'+'''All Yade displacements prescribed as boundary conditions
NonLinearStatic nsteps 2 contextOutputStep 1 rtolv 1.e-2 stiffMode 2 maxiter 50 controllmode 1 nmodules 0
domain 3dshell
OutputManager tstep_all dofman_all element_all
''')
inters=[i for i in O.interactions if (i.geom and i.phys)]
f.write("ndofman %d nelem %d ncrosssect 1 nmat 1 nbc %d nic 0 nltf 1 nbarrier 0\n"%(len(O.bodies),len(inters),len(O.bodies)*6))
bcMax=0; bcMap={}
for b in O.bodies:
mf=' '.join([str(a) for a in list(O.actions.f(b.id))+list(O.actions.m(b.id))])
f.write("## #%d: forces %s\n"%(b.id+1,mf))
f.write("Particle %d coords 3 %.15e %.15e %.15e rad %g"%(b.id+1,b.phys['se3'][0],b.phys['se3'][1],b.phys['se3'][2],b.mold['radius']))
bcMap[b.id]=tuple([bcMax+i for i in [1,2,3,4,5,6]])
bcMax+=6
f.write(' bc '+' '.join([str(i) for i in bcMap[b.id]])+'\n')
for j,i in enumerate(inters):
epsTNorm=sqrt(sum([e**2 for e in i.phys['epsT']]))
epsT='epsT [ %g %g %g ] %g'%(i.phys['epsT'][0],i.phys['epsT'][1],i.phys['epsT'][2],epsTNorm)
en=i.phys['epsN']; n=i.geom['normal']
epsN='epsN [ %g %g %g ] %g'%(en*n[0],en*n[1],en*n[2],en)
Fn='Fn [ %g %g %g ] %g'%(i.phys['normalForce'][0],i.phys['normalForce'][1],i.phys['normalForce'][2],i.phys['Fn'])
Fs='Fs [ %lf %lf %lf ] %lf'%(i.phys['shearForce'][0],i.phys['shearForce'][1],i.phys['shearForce'][2],sqrt(sum([fs**2 for fs in i.phys['shearForce']])))
f.write('## #%d #%d: %s %s %s %s\n'%(i.id1+1,i.id2+1,epsN,epsT,Fn,Fs))
f.write('CohSur3d %d nodes 2 %d %d mat 1 crossSect 1 area %g\n'%(j+1,i.id1+1,i.id2+1,i.phys['crossSection']))
# crosssection
f.write("SimpleCS 1 thick 1.0 width 1.0\n")
# material
ph=inters[0].phys
f.write("CohInt 1 kn %g ks %g e0 %g ef %g c 0. ksi %g coh %g tanphi %g d 1.0 conf 0.0 maxdist 0.0\n"%(ph['E'],ph['G'],ph['epsCrackOnset'],ph['epsFracture'],ph['xiShear'],ph['undamagedCohesion'],ph['tanFrictionAngle']))
# boundary conditions
for b in O.bodies:
displ=b.phys.displ; rot=b.phys.rot
dofs=[displ[0],displ[1],displ[2],rot[0],rot[1],rot[2]]
f.write('# particle %d\n'%b.id)
for dof in range(6):
f.write('BoundaryCondition %d loadTimeFunction 1 prescribedvalue %.15e\n'%(bcMap[b.id][dof],dofs[dof]))
#f.write('PiecewiseLinFunction 1 npoints 3 t 3 0. 10. 1000. f(t) 3 0. 10. 1000.\n')
f.write('ConstantFunction 1 f(t) 1.0\n')
[docs]def oofemDirectExport(fileBase,title=None,negIds=[],posIds=[]):
from yade.wrapper import Omega
material,bodies,interactions=[],[],[]
o=Omega()
strainers=[e for e in o.engines if e.name=='UniaxialStrainer']
if len(strainers)>0:
strainer=strainers[0]
posIds,negIds=strainer.posIds,strainer.negIds
else: strainer=None
f=open(fileBase+'.in','w')
# header
f.write(fileBase+'.out\n')
f.write((title if title else 'Yade simulation for '+fileBase)+'\n')
f.write("NonLinearStatic nsteps 2 contextOutputStep 1 rtolv 1.e-2 stiffMode 2 maxiter 50 controllmode 1 nmodules 0\n")
f.write("domain 3dShell\n")
f.write("OutputManager tstep_all dofman_all element_all\n")
inters=[i for i in o.interactions if (i.geom and i.phys)]
f.write("ndofman %d nelem %d ncrosssect 1 nmat 1 nbc 2 nic 0 nltf 1 nbarrier 0\n"%(len(o.bodies)+2,len(inters)))
# elements
f.write("Node 1 coords 3 0.0 0.0 0.0 bc 6 1 1 1 1 1 1\n")
f.write("Node 2 coords 3 0.0 0.0 0.0 bc 6 1 2 1 1 1 1\n")
for b in o.bodies:
f.write("Particle %d coords 3 %g %g %g rad %g"%(b.id+3,b.state.refPos[0],b.state.refPos[1],b.state.refPos[2],b.shape.radius))
if b.id in negIds: f.write(" dofType 6 1 1 1 1 1 1 masterMask 6 0 1 0 0 0 0 ")
elif b.id in posIds: f.write(" dofType 6 1 1 1 1 1 1 masterMask 6 0 2 0 0 0 0 0")
f.write('\n')
j=1
for i in inters:
f.write('CohSur3d %d nodes 2 %d %d mat 1 crossSect 1 area %g\n'%(j,i.id1+3,i.id2+3,i.phys.crossSection))
j+=1
# crosssection
f.write("SimpleCS 1 thick 1.0 width 1.0\n")
# material
ph=inters[0].phys
f.write("CohInt 1 kn %g ks %g e0 %g ef %g c 0. ksi %g coh %g tanphi %g damchartime %g damrateexp %g plchartime %g plrateexp %g d 1.0\n"%(ph.E,ph.G,ph.epsCrackOnset,ph.epsFracture,0.0,ph.undamagedCohesion,ph.tanFrictionAngle,ph.dmgTau,ph.dmgRateExp,ph.plTau,ph.plRateExp))
# boundary conditions
f.write('BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0\n')
f.write('BoundaryCondition 2 loadTimeFunction 1 prescribedvalue 1.e-4\n')
f.write('PiecewiseLinFunction 1 npoints 3 t 3 0. 10. 1000. f(t) 3 0. 10. 1000.\n')
[docs]def displacementsInteractionsExport(fName):
f=open(fName,'w')
print "Writing body displacements and interaction strains."
o=Omega()
for b in o.bodies:
x0,y0,z0=b.phys['refSe3'][0:3]; x,y,z=b.phys.pos
rx,ry,rz,rr=b.phys['se3'][3:7]
f.write('%d xyz [ %g %g %g ] dxyz [ %g %g %g ] rxyz [ %g %g %g ] \n'%(b.id,x0,y0,z0,x-x0,y-y0,z-z0,rr*rx,rr*ry,rr*rz))
f.write('\n')
for i in o.interactions:
if not i['isReal']:continue
epsTNorm=sqrt(sum([e**2 for e in i.phys['epsT']]))
epsT='epsT [ %g %g %g ] %g'%(i.phys['epsT'][0],i.phys['epsT'][1],i.phys['epsT'][2],epsTNorm)
en=i.phys['epsN']; n=i.geom['normal']
epsN='epsN [ %g %g %g ] %g'%(en*n[0],en*n[1],en*n[2],en)
Fn='Fn [ %g %g %g ] %g'%(i.phys['normalForce'][0],i.phys['normalForce'][1],i.phys['normalForce'][2],i.phys['Fn'])
Fs='Fs [ %lf %lf %lf ] %lf'%(i.phys['shearForce'][0],i.phys['shearForce'][1],i.phys['shearForce'][2],sqrt(sum([fs**2 for fs in i.phys['shearForce']])))
f.write('%d %d %s %s %s %s\n'%(i.id1,i.id2,epsN,epsT,Fn,Fs))
# f.write('%d %d %g %g %g %g %g\n'%(i.id1,i.id2,i.phys['epsN'],i.phys['epsT'][0],i.phys['epsT'][1],i.phys['epsT'][2],epsTNorm))
f.close()
[docs]def eliminateJumps(eps,sigma,numSteep=10,gapWidth=5,movWd=40):
from matplotlib.mlab import movavg
from numpy import diff,abs
import numpy
# get histogram of 'derivatives'
ds=abs(diff(sigma))
n,bins=numpy.histogram(ds)
i=1; sum=0
# numSteep steepest pieces will be discarded
while sum<numSteep:
#print n[-i],bins[-i]
sum+=n[-i]; i+=1
#print n[-i],bins[-i]
threshold=bins[-i]
# old algo: replace with nan's
##rEps,rSigma=eps[:],sigma[:]; nan=float('nan')
##indices=where(ds>threshold)[0]
##for i in indices:
## for ii in range(max(0,i-gapWidth),min(len(rEps),i+gapWidth+1)): rEps[ii]=rSigma[ii]=nan
## doesn't work with older numpy (?)
# indices1=where(ds>threshold)[0]
indices1=[]
for i in range(len(ds)):
if ds[i]>threshold: indices1.append(i)
indices=[]
for i in indices1:
for ii in range(i-gapWidth,i+gapWidth+1): indices.append(ii)
#print indices1, indices
rEps=[eps[i] for i in range(0,len(eps)) if i not in indices]
rSigma=[sigma[i] for i in range(0,len(sigma)) if i not in indices]
# apply moving average to the result
rSigma=movavg(rSigma,movWd)
rEps=rEps[movWd/2:-movWd/2+1]
return rEps,rSigma.flatten().tolist()