yade-dev team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

I wonder if this can be related to problems with O.load(), because
here I am NOT using O.load. I am only generating a new sample using
somewhat more complex method, including a call to O.reset().

I am attaching a stripped down script, but still a bit complex.

After I run that script, in the end I can get a following thing in python:

Yade [1]: O.bodies[0].state.dict()
 ->  [1]: 
{'accel': Vector3(0,0,0),
 'angAccel': Vector3(0,0,0),
 'angMom': Vector3(0,0,0),
 'angVel': Vector3(0,0,0),
 'blockedDOFs': 0,
 'inertia': Vector3(1.680768072756119e-10,1.680768072756119e-10,1.680768072756119e-10),
 'mass': 9.7475048955972981e-05,
 'refOri': Quaternion((1,0,0),0),
 'refPos': Vector3(0.046777726741579144,0.064726888517427708,0.060929861505333097),
 'se3': (Vector3(0.046777726741579144,0.064726888517427708,0.060929861505333097),
         Quaternion((1,0,0),0)),
 'vel': Vector3(0,0,0)}

Yade [2]: O.run()
Yade [3]: 48474 FATAL yade.ThreadRunner /home/janek/20-Programowanie/10-cpp/50-Yade/Code/HEAD/trunk/core/ThreadRunner.cpp:31 run: Exception occured: 
Body #0 has velocity==NaN!


Now - I have just checked that velocity is 0,0,0. And for some reason
it is NaN. why?

-- 
Janek Kozicki                               http://janek.kozicki.pl/  |

#!/usr/local/bin/yade-trunk
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-

# Script to perform a typical triaxial test
#-----------------------------------------------------
# simulation parameters
#-----------------------------------------------------
QUICK=True

size			= (0.1,0.1,0.1)
young			= 30e9 # 30 GPa  ## real
number_of_spheres	= 5000
local_poisson		= 0.3
frictionAngle		= 30
bending_ratio		= 0.5
moment_limit		= 0.1
std_dev_radius		= 0.5
avg_radius		= 0.0025 # r=0.0025m = 2.5mm ; d_50 = 0.05 m = 5 mm
max_compression		= 0.3
iso_compaction		= 50000
lateral_confinement	= 50000
damping			= 0.3
max_spheres_in_clump=6
min_spheres_in_clump=4
STAGE=0

if QUICK: #young=10000000 ## faster!
	young			= 15000000
	number_of_spheres	= 500
	max_compression		= 0.01

#-----------------------------------------------------
# materials
#-----------------------------------------------------
def set_materials(moment_limit_param):
	O.materials.append(CohFrictMat(	young=young,
					poisson=local_poisson,
					frictionAngle=radians(0),  # is set to zero for generating sample with low porosity. triax will set it later to good value, search for frictionAngle
					density=2600,
					isCohesive=False,
					isBroken=True,
					rollingStiffnessCoefficient=bending_ratio,
					momentLimitCoefficient=moment_limit_param,
					label='granular_material'))
	
	O.materials.append(CohFrictMat(	young=young,
					poisson=local_poisson,
					frictionAngle=radians(0),
					density=2600,
					isCohesive=False,
					isBroken=True,
					rollingStiffnessCoefficient=bending_ratio,
					momentLimitCoefficient=moment_limit_param,
					label='frictionless_material'))

#-----------------------------------------------------
# geometry
#-----------------------------------------------------
def set_geometry(id_spheres_radius_loc=None , clumps_loc=None):
	from yade import pack
	import itertools
	import random
	from numpy import arange
	from yade import utils
	def rnd():
		return random.random()

	if(clumps_loc == None):
		if(id_spheres_radius_loc == None):	
			# add the cloud of spheres
			sphere_cloud=pack.SpherePack()
			sphere_cloud.makeCloud(	minCorner=(0,0,0),
					maxCorner=size,
					rMean=avg_radius,
					rRelFuzz=std_dev_radius,
					num=number_of_spheres,
					periodic=False,
					porosity=-1)
			O.bodies.append([utils.sphere(	center,
						rad,
						material='granular_material')
						for center,rad in sphere_cloud])
		else:
			O.bodies.append([utils.sphere(	center,
						rad,
						material='granular_material')
						for center,rad in id_spheres_radius_loc.values()])
	else: # generate clumps
		used_spheres=[]
	#	print id_spheres_radius_loc
		for thisc in clumps_loc:
			col=(rnd(),rnd(),rnd())
			clump_spheres=[]
			for sphere_id in thisc:
				#O.bodies[sphere_id].shape.color=color
	#			print 'appending',sphere_id
				cent,rad = id_spheres_radius_loc[sphere_id]
				used_spheres.append(sphere_id)
	#			print 'center ',cent,'radius ',rad
				new_sphere=utils.sphere(center=cent,radius=rad,color=col,material='granular_material')
				clump_spheres.append(new_sphere)
			O.bodies.appendClumped(clump_spheres)
		for sphere_id in id_spheres_radius_loc:
	#		print 'checking sphere ', sphere_id
			if sphere_id not in used_spheres:
	#			print 'is not used, adding'
				cent,rad = id_spheres_radius_loc[sphere_id]
				col=(rnd(),rnd(),rnd())
				new_sphere=utils.sphere(center=cent,radius=rad,color=col,material='granular_material')
				O.bodies.append(new_sphere)
	#		else:
	#			print 'is used'

	# add the walls
	walls=yade.utils.aabbWalls(	
					#extrema=((0,0,0),size), # auto
					thickness=.01,
					oversizeFactor=1.2,
					material='frictionless_material')
	wallIds_local=O.bodies.append(walls)
	
	return wallIds_local

#-----------------------------------------------------
# list of engines
#-----------------------------------------------------
def set_engines(use_moments,wallIds):
	triax_local=TriaxialCompressionEngine(
		wall_left_id=wallIds[0],
		wall_right_id=wallIds[1],
		wall_bottom_id=wallIds[2],
		wall_top_id=wallIds[3],
		wall_back_id=wallIds[4],
		wall_front_id=wallIds[5],
		internalCompaction=0,
		autoCompressionActivation=1,
		sigmaIsoCompaction=iso_compaction,
		sigmaLateralConfinement=lateral_confinement,
		max_vel=1,
		StabilityCriterion=0.01, ######## 0.3 ??
		strainRate=0.1,
		frictionAngleDegree=frictionAngle,

		# also can use this, to change globally friction angle on all interactions
		# TriaxialCompressionEngine().setContactProperties(radians(25)) 

		fixedPorosity=1, # ??
		epsilonMax=max_compression,
		autoUnload=1,
		autoStopSimulation=1,
		noFiles=1,
		maxMultiplier=1.01,
		finalMaxMultiplier=1.001,
		radiusControlInterval=10,
		stiffnessUpdateInterval=10,
		fixedPoroCompaction=0
	)
	
	O.engines=[
		ForceResetter(),
		BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
		InsertionSortCollider(nBins=5,sweepLength=avg_radius*0.05),
		InteractionDispatchers(
			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
			
			#[Ip2_FrictMat_FrictMat_FrictPhys()],
			#[Law2_ScGeom_FrictPhys_Basic()]
			
			[Ip2_2xCohFrictMat_CohFrictPhys()],
			[Law2_ScGeom_CohFrictPhys_ElasticPlastic(always_use_moment_law=use_moments,momentRotationLaw=use_moments)]
		),
		GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=50),
		triax_local,
		NewtonIntegrator(damping=damping),
	]
	return triax_local

def save_spheres_data_and_find_clumps():
	import random
	id_spheres_radius={}
	ids=[]
	
	for b in O.bodies:
		#print b.id
		if b.shape.name=="Sphere":
			id_spheres_radius[b.id]=(b.state.se3[0],b.shape.radius)
			ids.append(b.id)
	
	interactions=[]
	used_spheres=[]
	
	for intr in O.interactions:
		a=intr.id1
		b=intr.id2
		if((a in ids) and (b in ids) and				# both are spheres
			(a != b) and						# aren't equal (rather impossible?)
			(intr.interactionGeometry!=None) and	
			(intr.interactionPhysics!=None) and	# existing interaction
			(intr.iterMadeReal>0)):
				interactions.append((a,b))
	
	random.shuffle(ids)
	
	clumps=[]
	
	for center_of_clump in ids:
		thisclump=[]
		thisclump.append(center_of_clump)
		used_spheres.append(center_of_clump)
		for a,b in interactions:
			if((a==center_of_clump) and (b not in used_spheres)):
				thisclump.append(b)
				used_spheres.append(b)
			if((b==center_of_clump) and (a not in used_spheres)):
				thisclump.append(a)
				used_spheres.append(a)
		if( len(thisclump) >= min_spheres_in_clump ):
			clumps.append(thisclump)
		else:
			for used in thisclump:
				used_spheres.remove(used)
	
	#for thisc in clumps:
	#	color=(rnd(),rnd(),rnd())
	#	for a in thisc:
	#		O.bodies[a].shape.color=color

	return id_spheres_radius,clumps

#-----------------------------------------------------
# run! moment law
STAGE = 1
# spheres, moment law
#-----------------------------------------------------
set_materials(moment_limit)
walls=set_geometry()
triax=set_engines(True,walls)
O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away

O.run(100,True)
O.dt=-1
O.run()

import time

time.sleep(10)

while triax.currentState!=3:
	time.sleep(1)

print "----- initial compression done -----"
print "--------- generating clumps --------"
id_spheres_radius_glob,clumps_glob = save_spheres_data_and_find_clumps()

O.pause()

#-----------------------------------------------------
# run again! - clumps, moment law
STAGE = 3
# clumps, moment law
#-----------------------------------------------------
O.reset()
set_materials(moment_limit)
walls=set_geometry(id_spheres_radius_glob,clumps_glob)
triax=set_engines(True,walls)

O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away

print(triax.currentState)

print 'Now type O.run() and see a NaN error'