yade-dev team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

I looked at his commit with command:

bzr diff -r before:2640..2640 | kompare -o -

but I don't see how it could help with unbalanced force.
Maybe I don't understand something here... and I modified
TriaxialStressController.cpp, while, maybe the solution is to
calculate unbalanced force somewhere else?

I don't know how TriaxialStressController.cpp is supposed to call
UnbalancedForceCallbacks.cpp. But it's a callback, and maybe some
callback mechanism has slipped my attention.

Currently that triax_clup.py is not working (the force never gets
balanced), unless the patch is applied.

best regards
Janek Kozicki

Bruno Chareyre said:     (by the date of Mon, 28 Mar 2011 12:55:04 +0200)

> Hi Janek,
> 
> Vaclav commited a fix for clumps unbalanced forces in r2640 (according
> to the commit log).
> Was it not working?
> I'll have a look.
> 
> Bruno
> 
> On 25/03/11 21:39, Janek Kozicki wrote:
> > Hi Bruno,
> >
> > you are very careful about TriaxialStressController.cpp, so I think
> > that first I will tell you what I changed, instead of surprising you with
> > a commit. So that you can check this yourself.
> >
> > I am attaching a small triax_clump.py which runs triaxial test on
> > clumps. This will not work, because unbalanced force never gets
> > balanced. Especially when you use "cigar" (line 82 in this script).
> >
> > The problem is described in detail in the attached patch, clump_unbalanced_force.diff
> > Apply this patch, and then the unbalanced force will handle clumps
> > correctly, I believe.
> >
> > If you like it, then commit this, else let me know what is wrong.
> >
> > best regards
> 
> 
> -- 
> _______________
> Bruno Chareyre
> Associate Professor
> ENSE³ - Grenoble INP
> Lab. 3SR
> BP 53 - 38041, Grenoble cedex 9 - France
> Tél : +33 4 56 52 86 21
> Fax : +33 4 76 82 70 43
> ________________
> 
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev
> More help   : https://help.launchpad.net/ListHelp
> 


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

=== modified file 'pkg/dem/TriaxialStressController.cpp'
--- pkg/dem/TriaxialStressController.cpp	2010-12-31 14:35:21 +0000
+++ pkg/dem/TriaxialStressController.cpp	2011-03-25 20:21:34 +0000
@@ -272,13 +272,41 @@
 		BodyContainer::iterator bi    = bodies->begin();
 		BodyContainer::iterator biEnd = bodies->end();
 		Real f;
+		typedef std::map<int,Vector3r> clump_forces_map;
+		clump_forces_map clump_forces;
 		for(  ; bi!=biEnd ; ++bi ) {
 			if ((*bi)->isDynamic()) {
-				f=getForce(scene,(*bi)->getId()).norm();
-				MeanUnbalanced += f;
-				if (f!=0) ++nBodies;
-			}
-		}
+				if((*bi)->isStandalone()) {
+					// this is correct only for stand alone spheres (non clump members)
+					f=getForce(scene,(*bi)->getId()).norm();
+					MeanUnbalanced += f;
+					if (f!=0) ++nBodies;
+				} else {
+					// for clumps I cannot add the norm() of forces of all clump members
+					// first I must sum all forces from all spheres in a clump, then take norm() of it
+					// because there are forces in a clump that do not cross through its centroid,
+					// but can get cancelled overall
+					int cid((*bi)->clumpId);
+					if(clump_forces.find(cid) == clump_forces.end()) {
+						// that's the first member of this clump
+						clump_forces[cid] = getForce(scene,(*bi)->getId());
+					} else { // there's already some force stored here, so add it
+						clump_forces[cid] += getForce(scene,(*bi)->getId());
+					}
+				}
+			}
+		}
+
+		if( !clump_forces.empty()) {
+			FOREACH(const clump_forces_map::value_type& v, clump_forces) {
+				f = v.second.norm();
+				if(f != 0) {
+					MeanUnbalanced += f;
+					++nBodies;
+				}
+			}
+		}
+
 		if (nBodies != 0 && MeanForce != 0) MeanUnbalanced = MeanUnbalanced/nBodies/MeanForce;
 		return  MeanUnbalanced;
 	} else {

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

# Script to perform a typical triaxial test
#-----------------------------------------------------
# simulation parameters
#-----------------------------------------------------
size			= Vector3(0.1,0.1,0.1)
young			= 3e8
number_of_spheres	= 1000
local_poisson		= 0.3
frictionAngle		= 10
bending_ratio		= 0.03
moment_limit		= 0.5
std_dev_radius		= 0.5
avg_radius		= 0.0025
max_compression		= 0.12
iso_compaction		= 200000
lateral_confinement	= 200000
damping			= 0.3
clump_DT_coeff		= 0.05
gen_clumps_grid_size	= 10
gen_clumps_grid_size	= (number_of_spheres**(0.33))*0.5

#-----------------------------------------------------
# materials
#-----------------------------------------------------
def set_materials(moment_limit_param,use_moments):
	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,
					alphaKr  =bending_ratio,
					alphaKtw =bending_ratio,
					momentRotationLaw=use_moments,
					etaRoll=moment_limit_param,
					label='granular_material'))
	
	O.materials.append(CohFrictMat(	young=young,
					poisson=local_poisson,
					frictionAngle=radians(0),
					density=2600,
					isCohesive=False,
					alphaKr  =bending_ratio,
					alphaKtw =bending_ratio,
					momentRotationLaw=use_moments,
					etaRoll=moment_limit_param,	
					label='frictionless_material'))

#-----------------------------------------------------
# geometry
#-----------------------------------------------------
def ball2(pos,scale):
	col=utils.randomColor()
	cigar_spheres=[]
	cigar_spheres.append(utils.sphere(center=pos+0.25*Vector3(0,0,scale),radius=0.5 *scale,color=col,material='granular_material'))
	cigar_spheres.append(utils.sphere(center=pos-0.25*Vector3(0,0,scale),radius=0.5 *scale,color=col,material='granular_material'))
	return cigar_spheres

def cigar(pos,scale):
	col=utils.randomColor()
	cigar_spheres=[]
	cigar_spheres.append(utils.sphere(center=pos                        ,radius=0.5 *scale,color=col,material='granular_material'))
	cigar_spheres.append(utils.sphere(center=pos+0.25*Vector3(0,0,scale),radius=0.48*scale,color=col,material='granular_material'))
	cigar_spheres.append(utils.sphere(center=pos-0.25*Vector3(0,0,scale),radius=0.48*scale,color=col,material='granular_material'))
	cigar_spheres.append(utils.sphere(center=pos+0.5 *Vector3(0,0,scale),radius=0.41*scale,color=col,material='granular_material'))
	cigar_spheres.append(utils.sphere(center=pos-0.5 *Vector3(0,0,scale),radius=0.41*scale,color=col,material='granular_material'))
	cigar_spheres.append(utils.sphere(center=pos+0.75*Vector3(0,0,scale),radius=0.25*scale,color=col,material='granular_material'))
	cigar_spheres.append(utils.sphere(center=pos-0.75*Vector3(0,0,scale),radius=0.25*scale,color=col,material='granular_material'))
	return cigar_spheres

def clump_grid(dx,dy,dz,d,mn,mx):
	import random
	from yade import pack
	import itertools
	from numpy import arange
	for xyz in itertools.product(arange(0,d),arange(0,d),arange(0,d)):
		ids_spheres =	O.bodies.appendClumped(
			#ball2(
			cigar(
				Vector3(mn[0]+xyz[0]*dx,mn[0]+xyz[1]*dy,mn[2]+xyz[2]*dz),
				avg_radius + std_dev_radius*avg_radius*(random.random()-0.5)
			))

def set_geometry(clumps_grid_size):
	from yade import utils
	# generate clumps
	print '================== clumps'
	d=clumps_grid_size
	mn,mx=Vector3(0.003,0.003,0.003),size
	dx=(mx[0]-mn[0])/d
	dy=(mx[1]-mn[1])/d
	dz=(mx[2]-mn[2])/d
	clump_grid(dx,dy,dz,d,mn,mx)
	clump_grid(dx,dy,dz,d,mn+(0.5*Vector3(dx,dy,dz)),mx+(0.5*Vector3(dx,dy,dz)))
	# add the walls
	walls=yade.utils.aabbWalls(	
					#extrema=((0,0,0),size), # auto
					thickness=.01,
					oversizeFactor=2.0,
					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.3,
		strainRate=0.1,
		frictionAngleDegree=frictionAngle,
		fixedPoroCompaction=False,
		#fixedPorosity=0.347, # ??
		fixedPorosity=0.4, # ??
		epsilonMax=max_compression,
		autoUnload=1,
		autoStopSimulation=1,
		noFiles=1,
		maxMultiplier=1.01,
		finalMaxMultiplier=1.001,
		radiusControlInterval=10,
		stiffnessUpdateInterval=10,
	)
	
	O.engines=[
		ForceResetter(),
		BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
		InsertionSortCollider(nBins=5,sweepLength=avg_radius*0.05),
		InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
			[Ip2_2xCohFrictMat_CohFrictPhys()],
			[Law2_ScGeom_CohFrictPhys_CohesionMoment(always_use_moment_law=use_moments)]
		),
		GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=50),
		triax_local,
		NewtonIntegrator(damping=damping),
		PyRunner(iterPeriod=1000,command='plot1()'),
	]
	return triax_local

#-----------------------------------------------------
# plot some results
#-----------------------------------------------------
from yade import plot

plot.plots={'epsilon_1':('sigma_1',None,'porosity','void_ratio'),'epsilon_1_':('epsilon_123')}

def plot1():
	if(triax.currentState==3):
		p=utils.voxelPorosityTriaxial(triax,300,avg_radius)
		plot.addData(
			sigma_1=triax.stress(1)[1],
			epsilon_1=triax.strain[1],
			epsilon_1_=triax.strain[1],
			epsilon_1__=triax.strain[1],
			epsilon_123=(-1.0)*(triax.volumetricStrain),
			porosity=p,
			void_ratio=p/(1-p)
		     )

#-----------------------------------------------------
# set gui
#-----------------------------------------------------
from yade import qt
qt.View()
qt.Controller()
O.reset()
set_materials(0,False)
walls=set_geometry(gen_clumps_grid_size)
triax=set_engines(False,walls)
O.dt=0.5*utils.PWaveTimeStep() # initial timestep, to not explode right away
yade.qt.controller.setWindowTitle('CLUMP TEST')
plot.plot()