yade-users team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

Hi

It's a fun coincidence that currently I'm adding aggregates to my
lattice model. So that I will simulate concrete with aggregates.

I just want to send you now my filegenerator for lattice, the code that
generates aggregates is at the bottom of LatticeExample.cpp
I use geometry-based generation method. Pick random coordinates and
random radius and try to put it into the speciemen. Do this until
certian aggregates compacticy is obtaioned, very primitive method. It is 2D.

Just compile this filegenerator, run it, then load scene.xml. Make sure
that you set SimulationController->Display->drawWireFrame option. Now
you can see generated aggragates (nodes that belong to an aggregate have
black color)

look into the code (bottom of file). You will see that I use
boost::random library - I recommend using it:

// random
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/normal_distribution.hpp>


// this below creates a random number generator between 0 and 1 with
// uniform distribution

static boost::variate_generator<StdGenerator&, boost::uniform_real<> >
		random1(generator, boost::uniform_real<>(0,1));


// this creates a random number generator with Gaussian distribution
// the parameter is mean, and sigma (parameters of gaussian distribution)

	static boost::variate_generator<StdGenerator&, boost::normal_distribution<> >
		randomN(generator, boost::normal_distribution<>(aggregateMeanRadius,aggregateSigmaRadius));


boost::random has more random distributions to choose from, look here:

http://www.boost.org/libs/random/index.html

http://www.boost.org/libs/random/random-distributions.html


I hope that this will help you with problem "I need to find a random number generator" :)
I you have any questions just ask.


-- 
Janek Kozicki                                                         |

/*************************************************************************
*  Copyright (C) 2004 by Janek Kozicki                                   *
*  cosurgi@xxxxxxxxxx                                                    *
*                                                                        *
*  This program is free software; it is licensed under the terms of the  *
*  GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/

#include "LatticeExample.hpp"

#include "LatticeSetParameters.hpp"
#include "LatticeBeamParameters.hpp"
#include "LatticeBeamAngularSpring.hpp"
#include "NonLocalDependency.hpp"
#include "LatticeNodeParameters.hpp"
#include "LineSegment.hpp"
#include "LatticeLaw.hpp"
//#include "StrainRecorder.hpp"
//#include "NonLocalInitializer.hpp"

#include <yade/yade-package-common/Sphere.hpp>

#include <yade/yade-core/Body.hpp>
#include <yade/yade-package-common/MetaInteractingGeometry.hpp>
#include <yade/yade-package-common/BoundingVolumeMetaEngine.hpp>
#include <yade/yade-package-common/GeometricalModelMetaEngine.hpp>

#include <yade/yade-package-common/AABB.hpp>

#include <yade/yade-package-common/BodyRedirectionVector.hpp>
#include <yade/yade-package-common/InteractionVecSet.hpp>
#include <yade/yade-package-common/PhysicalActionVectorVector.hpp>

#include <yade/yade-package-common/DisplacementEngine.hpp>
#include <yade/yade-package-common/PhysicalParametersMetaEngine.hpp>
#include <yade/yade-package-common/PhysicalActionApplier.hpp>

#include <yade/yade-package-common/PhysicalActionContainerInitializer.hpp>


using namespace boost;
using namespace std;


LatticeExample::LatticeExample() : FileGenerator()
{
	nodeGroupMask 		= 1;
	beamGroupMask 		= 2;
	
	speciemen_size_in_meters = Vector3r(0.1,0.1,0.0001);
	cellsizeUnit_in_meters	 = 0.002;
	minAngle_betweenBeams_deg= 20.0;
	disorder_in_cellsizeUnit = Vector3r(0.6,0.6,0.0);
	maxLength_in_cellsizeUnit= 1.9;
	triangularBaseGrid 	 = true;
	useNonLocalModel 	 = false;
	useBendTensileSoftening	 = false;
	useStiffnessSoftening	 = false;
	nonLocalL_in_cellsizeUnit = 1.1; 	// l
				
	crit_TensileStrain_percent = 100.0;	// E_min
	crit_ComprStrain_percent   = 50.0;	// E_max
			
	longitudalStiffness_noUnit = 1.0;	// k_l
	bendingStiffness_noUnit    = 0.0;	// k_b
	
	region_A_min 		= Vector3r(-0.006, 0.096,-1);
	region_A_max 		= Vector3r( 0.16 , 0.16 , 1);
	direction_A 		= Vector3r(0,1,0);
	displacement_A_meters	= 0.0001;
	
	region_B_min 		= Vector3r(-0.006,-0.006,-1);
	region_B_max 		= Vector3r( 0.16 , 0.004, 1);
	direction_B 		= Vector3r(0,-1,0);
	displacement_B_meters	= 0.0001;

	region_C_min 		= Vector3r(-0.006, 0.096,-1);
	region_C_max 		= Vector3r( 0.16 , 0.16 , 1);
	direction_C 		= Vector3r(0,1,0);
	displacement_C_meters	= 0.0001;
	
	region_D_min 		= Vector3r(-0.006,-0.006,-1);
	region_D_max 		= Vector3r( 0.16 , 0.004, 1);
	direction_D 		= Vector3r(0,-1,0);
	displacement_D_meters	= 0.0001;
	
	strainRecorder_xz_plane	= -1;
	outputFile 		= "../data/strains";
	subscribedBodies.clear();
		
	regionDelete_A_min 	= Vector3r(0,0,0);
	regionDelete_A_max 	= Vector3r(0,0,0);
	regionDelete_B_min 	= Vector3r(0,0,0);
	regionDelete_B_max 	= Vector3r(0,0,0);
	regionDelete_C_min 	= Vector3r(0,0,0);
	regionDelete_C_max 	= Vector3r(0,0,0);
	regionDelete_D_min 	= Vector3r(0,0,0);
	regionDelete_D_max 	= Vector3r(0,0,0);
	regionDelete_E_min 	= Vector3r(0,0,0);
	regionDelete_E_max 	= Vector3r(0,0,0);
	regionDelete_F_min 	= Vector3r(0,0,0);
	regionDelete_F_max 	= Vector3r(0,0,0);

	nonDestroy_A_min	= Vector3r(0,0,0);
	nonDestroy_A_max	= Vector3r(0,0,0);
	nonDestroy_B_min	= Vector3r(0,0,0);
	nonDestroy_B_max	= Vector3r(0,0,0);

	useAggregates		= true;
	aggregatePercent	= 60;
	aggregateMeanRadius	= cellsizeUnit_in_meters*3;
	aggregateSigmaRadius	= cellsizeUnit_in_meters;
	agg_longStiffness_noUnit= longitudalStiffness_noUnit*10;	// k_l
	agg_bendStiffness_noUnit= bendingStiffness_noUnit*10;		// k_b
}


LatticeExample::~LatticeExample()
{

}


void LatticeExample::registerAttributes()
{
	REGISTER_ATTRIBUTE(speciemen_size_in_meters); 	// size
	REGISTER_ATTRIBUTE(cellsizeUnit_in_meters);	// g [m]  	- cell size
	REGISTER_ATTRIBUTE(minAngle_betweenBeams_deg); 	// a [deg] 	- min angle
	REGISTER_ATTRIBUTE(disorder_in_cellsizeUnit); 	// s [-] 	- disorder 
	REGISTER_ATTRIBUTE(maxLength_in_cellsizeUnit);	// r [-] 	- max beam length
	
	REGISTER_ATTRIBUTE(crit_TensileStrain_percent); // E_min [%]    - default 0.02 %
	REGISTER_ATTRIBUTE(crit_ComprStrain_percent);   // E_max [%]    - default 0.2 %
	REGISTER_ATTRIBUTE(longitudalStiffness_noUnit); // k_l [-]      - default 1.0
	REGISTER_ATTRIBUTE(bendingStiffness_noUnit);    // k_b [-]      - default 0.6
	
	REGISTER_ATTRIBUTE(triangularBaseGrid); 	// 		- triangles
	REGISTER_ATTRIBUTE(useBendTensileSoftening);
	REGISTER_ATTRIBUTE(useStiffnessSoftening);
	REGISTER_ATTRIBUTE(useNonLocalModel);
	REGISTER_ATTRIBUTE(nonLocalL_in_cellsizeUnit); 	// l
	
	REGISTER_ATTRIBUTE(region_A_min);
	REGISTER_ATTRIBUTE(region_A_max);
	REGISTER_ATTRIBUTE(direction_A);
	REGISTER_ATTRIBUTE(displacement_A_meters);
	
	REGISTER_ATTRIBUTE(region_B_min);
	REGISTER_ATTRIBUTE(region_B_max);
	REGISTER_ATTRIBUTE(direction_B);
	REGISTER_ATTRIBUTE(displacement_B_meters);
	
	REGISTER_ATTRIBUTE(region_C_min);
	REGISTER_ATTRIBUTE(region_C_max);
	REGISTER_ATTRIBUTE(direction_C);
	REGISTER_ATTRIBUTE(displacement_C_meters);
	
	REGISTER_ATTRIBUTE(region_D_min);
	REGISTER_ATTRIBUTE(region_D_max);
	REGISTER_ATTRIBUTE(direction_D);
	REGISTER_ATTRIBUTE(displacement_D_meters);
	
	REGISTER_ATTRIBUTE(strainRecorder_xz_plane);
	REGISTER_ATTRIBUTE(outputFile);
	
	REGISTER_ATTRIBUTE(regionDelete_A_min);
	REGISTER_ATTRIBUTE(regionDelete_A_max);
	REGISTER_ATTRIBUTE(regionDelete_B_min);
	REGISTER_ATTRIBUTE(regionDelete_B_max);
	REGISTER_ATTRIBUTE(regionDelete_C_min);
	REGISTER_ATTRIBUTE(regionDelete_C_max);
	REGISTER_ATTRIBUTE(regionDelete_D_min);
	REGISTER_ATTRIBUTE(regionDelete_D_max);
	REGISTER_ATTRIBUTE(regionDelete_E_min);
	REGISTER_ATTRIBUTE(regionDelete_E_max);
	REGISTER_ATTRIBUTE(regionDelete_F_min);
	REGISTER_ATTRIBUTE(regionDelete_F_max);

	REGISTER_ATTRIBUTE(nonDestroy_A_min);
	REGISTER_ATTRIBUTE(nonDestroy_A_max);
	REGISTER_ATTRIBUTE(nonDestroy_B_min);
	REGISTER_ATTRIBUTE(nonDestroy_B_max);

	REGISTER_ATTRIBUTE(useAggregates);
	REGISTER_ATTRIBUTE(aggregatePercent);
	REGISTER_ATTRIBUTE(aggregateMeanRadius);
	REGISTER_ATTRIBUTE(aggregateSigmaRadius);
	REGISTER_ATTRIBUTE(agg_longStiffness_noUnit);
	REGISTER_ATTRIBUTE(agg_bendStiffness_noUnit);
}

string LatticeExample::generate()
{
	rootBody = shared_ptr<MetaBody>(new MetaBody);
	createActors(rootBody);
	positionRootBody(rootBody);

	
	rootBody->persistentInteractions	= shared_ptr<InteractionContainer>(new InteractionVecSet);
	rootBody->volatileInteractions		= shared_ptr<InteractionContainer>(new InteractionVecSet);
	rootBody->actionParameters		= shared_ptr<PhysicalActionContainer>(new PhysicalActionVectorVector);
	rootBody->bodies 			= shared_ptr<BodyContainer>(new BodyRedirectionVector);

	
	shared_ptr<Body> body;
	
	Vector3r nbNodes = speciemen_size_in_meters / cellsizeUnit_in_meters;
	if(triangularBaseGrid)
		nbNodes[1] *= 1.15471; // bigger by sqrt(3)/2 factor

	unsigned int totalNodesCount = 0;

	for( int j=0 ; j<=nbNodes[1] ; j++ )
	{
		for( int i=0 ; i<=nbNodes[0] ; i++ )
			for( int k=0 ; k<=nbNodes[2] ; k++)
			{
				shared_ptr<Body> node;
				if(createNode(node,i,j,k))
					rootBody->bodies->insert(node), ++totalNodesCount;
			}
			
		// strain recorder - set what to monitor
		if( j == 0 && subscribedBodies.size() == 0 )
			subscribedBodies.push_back(static_cast<unsigned int>(nbNodes[0]) / 2);
	//	if( j == static_cast<int>(nbNodes[0]) && subscribedBodies.size() == 1 )
	//		subscribedBodies.push_back(totalNodesCount- 4);// - static_cast<int>(nbNodes[0]) / 2);
		//	subscribedBodies.push_back(  	  (static_cast<int>(nbNodes[0]))
		//					* (static_cast<int>(nbNodes[1]))
		//					* (static_cast<int>(nbNodes[2]))
		//					- static_cast<int>(nbNodes[0])/2);
	}
	subscribedBodies.push_back(totalNodesCount - static_cast<int>(nbNodes[0]) / 2);
	
	BodyRedirectionVector bc;
	bc.clear();

	BodyContainer::iterator bi    = rootBody->bodies->begin();
	BodyContainer::iterator bi2;
	BodyContainer::iterator biEnd = rootBody->bodies->end();
	int beam_counter = 0;
	float nodes_a=0;
	float nodes_all = rootBody->bodies->size();
	for(  ; bi!=biEnd ; ++bi )  // loop over all nodes, to create beams
	{
		Body* bodyA = (*bi).get(); // first_node
	
		bi2 = bi;
		++bi2;
		nodes_a+=1.0;
		
		for( ; bi2!=biEnd ; ++bi2 )
		{
			Body* bodyB = (*bi2).get(); // all other nodes
			// warning - I'm assuming that there are ONLY Nodes in the rootBody
			LatticeNodeParameters* a = static_cast<LatticeNodeParameters*>(bodyA->physicalParameters.get());
			LatticeNodeParameters* b = static_cast<LatticeNodeParameters*>(bodyB->physicalParameters.get());
			
			if ((a->se3.position - b->se3.position).squaredLength() < std::pow(maxLength_in_cellsizeUnit*cellsizeUnit_in_meters,2) )  
			{
				shared_ptr<Body> beam;
				createBeam(beam,bodyA->getId(),bodyB->getId());
				calcBeamPositionOrientationLength(beam);
				if(checkMinimumAngle(bc,beam))
				{
					if( ++beam_counter % 100 == 0 )
						cerr << "creating beam: " << beam_counter << " , " << ((nodes_a/nodes_all)*100.0)  << " %\n"; 
					
					bc.insert(beam);
				}
			}
		}
	}

	bi    = bc.begin();
	biEnd = bc.end();
	for(  ; bi!=biEnd ; ++bi )  // loop over all newly created beams ...
	{
		shared_ptr<Body> b = *bi;
		rootBody->bodies->insert(b); // .. to insert then into rootBody
		
		// attach them to strain recorder
		
		Real xz_plane = strainRecorder_xz_plane;
		Real pos1 = (*(rootBody->bodies))[static_cast<LatticeBeamParameters*>( b->physicalParameters.get() )->id1]->physicalParameters->se3.position[1]; // beam first node
		Real pos2 = (*(rootBody->bodies))[static_cast<LatticeBeamParameters*>( b->physicalParameters.get() )->id2]->physicalParameters->se3.position[1]; // beam second node
		if( pos1 < pos2 )
			std::swap(pos1,pos2); // make sure that pos1 is bigger 
		if(        pos1 > xz_plane
			&& pos2 < xz_plane 
			)
			subscribedBodies.push_back( b->getId() ); // beam crosses the plane!
	}
	//strainRecorder->subscribedBodies 	= subscribedBodies;
	//strainRecorder->initialLength =   (*(rootBody->bodies))[subscribedBodies[0]]->physicalParameters->se3.position[1] 
	//				- (*(rootBody->bodies))[subscribedBodies[1]]->physicalParameters->se3.position[1];
		

	{ // remember what node is in contact with what beams
		//    node                   beams
		connections.resize(totalNodesCount);
		
		bi    = rootBody->bodies->begin();
		biEnd = rootBody->bodies->end();
		
		for(  ; bi!=biEnd ; ++bi )  // loop over all beams
		{
			Body* body = (*bi).get();
			if( ! ( body->getGroupMask() & beamGroupMask ) )
				continue; // skip non-beams
			
			LatticeBeamParameters* beam = static_cast<LatticeBeamParameters*>(body->physicalParameters.get() );
			connections[beam->id1].push_back(body->getId());
			connections[beam->id2].push_back(body->getId());
		}
	}
	
	{ // create angular springs between beams
		bi    = rootBody->bodies->begin();
		biEnd = rootBody->bodies->end();
		float all_bodies = rootBody->bodies->size();
		int current = 0;
		for(  ; bi!=biEnd ; ++bi )  // loop over all beams
		{
			if( ++current % 100 == 0 )
				cerr << "angular springs: " << current << " , " << ((static_cast<float>(current)/all_bodies)*100.0) << " %\n";
				
			Body* body = (*bi).get();
			if( ! ( body->getGroupMask() & beamGroupMask ) )
				continue; // skip non-beams
				
			calcBeamAngles(body,rootBody->bodies.get(),rootBody->persistentInteractions.get());
		}
	}
	
	regionDelete(rootBody,regionDelete_A_min,regionDelete_A_max);
	regionDelete(rootBody,regionDelete_B_min,regionDelete_B_max);
	regionDelete(rootBody,regionDelete_C_min,regionDelete_C_max);
	regionDelete(rootBody,regionDelete_D_min,regionDelete_D_max);
	regionDelete(rootBody,regionDelete_E_min,regionDelete_E_max);
	regionDelete(rootBody,regionDelete_F_min,regionDelete_F_max);

	nonDestroy(rootBody,nonDestroy_A_min,nonDestroy_A_max);
	nonDestroy(rootBody,nonDestroy_B_min,nonDestroy_B_max);

	imposeTranslation(rootBody,region_A_min,region_A_max,direction_A,displacement_A_meters);
	imposeTranslation(rootBody,region_B_min,region_B_max,direction_B,displacement_B_meters);
	imposeTranslation(rootBody,region_C_min,region_C_max,direction_C,displacement_C_meters);
	imposeTranslation(rootBody,region_D_min,region_D_max,direction_D,displacement_D_meters);

	if(useAggregates) addAggregates(rootBody);
	
	cerr << "finished.. saving\n";

 	return "Number of nodes created:\n" + lexical_cast<string>(nbNodes[0]) + ","
	 				    + lexical_cast<string>(nbNodes[1]) + ","
					    + lexical_cast<string>(nbNodes[2]);

}

/// returns true if angle is bigger than minAngle_betweenBeams_deg
bool LatticeExample::checkAngle( Vector3r a, Vector3r& b)
{
	Quaternionr al;
	al.align(a,b);
	Vector3r axis;
	Real angle;
	al.toAxisAngle(axis, angle);
	angle *= 180.0/Mathr::PI ;
//	cerr << " angle: " << angle << "\n";
	return angle > minAngle_betweenBeams_deg;
}

/// returns true if angle is bigger than minAngle_betweenBeams_deg
bool LatticeExample::checkMinimumAngle(BodyRedirectionVector& bc,shared_ptr<Body>& body)
{
	bool answer = true;
	LatticeBeamParameters* newBeam = static_cast<LatticeBeamParameters*>(body->physicalParameters.get());
	
	BodyContainer::iterator bi    = bc.begin();
	BodyContainer::iterator biEnd = bc.end();
	for(  ; bi!=biEnd ; ++bi )  // loop over all beams - they MUST be beams, for static_cast<> 
	{ 
		LatticeBeamParameters* oldBeam = static_cast<LatticeBeamParameters*>((*bi)->physicalParameters.get());
		if( 	   (oldBeam->id1 == newBeam->id1)
			|| (oldBeam->id2 == newBeam->id2))
			answer = answer && checkAngle(   oldBeam->direction ,  newBeam->direction );
		if( 	   (oldBeam->id2 == newBeam->id1)
			|| (oldBeam->id1 == newBeam->id2))
			answer = answer && checkAngle( - oldBeam->direction ,  newBeam->direction );
	} 
	return answer;
}

bool LatticeExample::createNode(shared_ptr<Body>& body, int i, int j, int k)
{
	body = shared_ptr<Body>(new Body(0,nodeGroupMask));
	shared_ptr<LatticeNodeParameters> physics(new LatticeNodeParameters);
	shared_ptr<Sphere> gSphere(new Sphere);
	
	Quaternionr q;
	q.fromAxisAngle( Vector3r(Mathr::unitRandom(),Mathr::unitRandom(),Mathr::unitRandom()) , Mathr::unitRandom()*Mathr::PI );
	
	float  triang_x = triangularBaseGrid ? (static_cast<float>(j%2))*0.5 : 0;
	double triang_y = triangularBaseGrid ? 0.86602540378443864676        : 1; // sqrt(3)/2
	
	Vector3r position		= ( Vector3r(i+triang_x,j*triang_y,k)
					  + Vector3r( 	  Mathr::symmetricRandom()*disorder_in_cellsizeUnit[0]
					  		, Mathr::symmetricRandom()*disorder_in_cellsizeUnit[1]
							, Mathr::symmetricRandom()*disorder_in_cellsizeUnit[2]
						    ) * 0.5 // *0.5 because symmetricRandom is (-1,1), and disorder is whole size where nodes can appear
					  )*cellsizeUnit_in_meters;

	if( 	   position[0] >= speciemen_size_in_meters[0] 
		|| position[1] >= speciemen_size_in_meters[1]
		|| position[2] >= speciemen_size_in_meters[2] )
		return false;

	Real radius 			= cellsizeUnit_in_meters*0.05;
	
	body->isDynamic			= true;
	
	physics->se3			= Se3r(position,q);

	gSphere->radius			= radius;
	gSphere->diffuseColor		= Vector3f(0.8,0.8,0.8);
	gSphere->wire			= false;
	gSphere->visible		= true;
	gSphere->shadowCaster		= false;
	
	body->geometricalModel		= gSphere;
	body->physicalParameters	= physics;

//////////////////////
// do tego przyk³adu ustawiam wspó³rzêdne - FIXME - przyda³oby siê ³adowanie wspó³rzêdnych wierzcho³ków sk±din±d
//static int zzzzz=0;
//switch(zzzzz%5)
//{
//	case 0 	: physics->se3.position=Vector3r(0.4,1.5,0); break;
//	case 1 	: physics->se3.position=Vector3r(0.8,0.6,0); break;
//	case 2 	: physics->se3.position=Vector3r(0,0,0); break;
//	case 3 	: physics->se3.position=Vector3r(1.6,0.5,0); break;
//	case 4 	: physics->se3.position=Vector3r(2.0,0,0); break;
//};
//zzzzz++;
//
//////////////////////
	
	return true;
}


void LatticeExample::createBeam(shared_ptr<Body>& body, unsigned int i, unsigned int j)
{
	body = shared_ptr<Body>(new Body(0,beamGroupMask));
	shared_ptr<LatticeBeamParameters> physics(new LatticeBeamParameters);
	shared_ptr<LineSegment> gBeam(new LineSegment);
	
	Real length 			= 1.0; // unspecified for now, calcBeamsPositionOrientationLength will calculate it
	
	body->isDynamic			= true;
	
	physics->id1 			= i;
	physics->id2 			= j;

	gBeam->length			= length;
	gBeam->diffuseColor		= Vector3f(0.6,0.6,0.6);
	gBeam->wire			= false;
	gBeam->visible			= true;
	gBeam->shadowCaster		= false;
	
	body->geometricalModel		= gBeam;
	body->physicalParameters	= physics;
}


void LatticeExample::calcBeamPositionOrientationLength(shared_ptr<Body>& body)
{
	LatticeBeamParameters* beam = static_cast<LatticeBeamParameters*>(body->physicalParameters.get());
	shared_ptr<Body>& bodyA = (*(rootBody->bodies))[beam->id1];
	shared_ptr<Body>& bodyB = (*(rootBody->bodies))[beam->id2];
	Se3r& se3A 		= bodyA->physicalParameters->se3;
	Se3r& se3B 		= bodyB->physicalParameters->se3;
	
	Se3r se3Beam;
	se3Beam.position 	= (se3A.position + se3B.position)*0.5;
	Vector3r dist 		= se3A.position - se3B.position;
	
	Real length 		= dist.normalize();
	beam->direction 	= dist;
	beam->initialDirection 	= dist;
	beam->length 		= length;
	beam->initialLength 	= length;
	
	beam->criticalTensileStrain 	= crit_TensileStrain_percent/100.0;
	beam->criticalCompressiveStrain = crit_ComprStrain_percent/100.0;
	beam->longitudalStiffness 	= longitudalStiffness_noUnit;
	beam->bendingStiffness 		= bendingStiffness_noUnit;
	
	se3Beam.orientation.align( Vector3r::UNIT_X , dist );
	beam->se3 		= se3Beam;
	beam->se3Displacement.position 	= Vector3r(0.0,0.0,0.0);
	beam->se3Displacement.orientation.align(dist,dist);
}

void LatticeExample::calcAxisAngle(LatticeBeamParameters* beam, BodyContainer* bodies, unsigned int otherId, InteractionContainer* ints, unsigned int thisId)
{ 
	if( ! ints->find(otherId,thisId) && otherId != thisId )
	{
		LatticeBeamParameters* 	otherBeam 		= static_cast<LatticeBeamParameters*>( ((*(bodies))[ otherId ])->physicalParameters.get() );
		
	//	Quaternionr 		rotation;
	//	Vector3r 		axis;
		Real 			angle;
		
	//	rotation.align( beam->direction , otherBeam->direction );
	//	rotation.toAxisAngle (axis, angle);
	//	Vector3r result = axis*angle;
	
		angle = (beam->direction.cross(otherBeam->direction)[2] > 0.0 ? 1.0 : -1.0) * Mathr::aCos( beam->direction.dot(otherBeam->direction) );

		shared_ptr<Interaction> 		interaction(new Interaction( thisId , otherId ));
		shared_ptr<LatticeBeamAngularSpring> 	angularSpring(new LatticeBeamAngularSpring);
		
	//	angularSpring->initialAngle 		= result;
	//	angularSpring->angle 			= result;
	//	angularSpring->initialAngle 		= rotation;
	//	angularSpring->angle 			= rotation;
		angularSpring->initialAngle 		= angle;
		angularSpring->angle 			= angle;
		
		interaction->isReal			= true;
		interaction->isNew 			= false;
		interaction->interactionPhysics 	= angularSpring;
		ints->insert(interaction);
	}
}

void LatticeExample::calcBeamAngles(Body* body, BodyContainer* bodies, InteractionContainer* ints)
{
	LatticeBeamParameters* beam 	= static_cast<LatticeBeamParameters*>(body->physicalParameters.get());

	std::vector<unsigned int>::iterator i   = connections[beam->id1].begin();
	std::vector<unsigned int>::iterator end = connections[beam->id1].end();
	
	for( ; i != end ; ++i )
		calcAxisAngle(beam,bodies,*i,ints,body->getId());
	
	i   = connections[beam->id2].begin();
	end = connections[beam->id2].end();
	for( ; i != end ; ++i )
		calcAxisAngle(beam,bodies,*i,ints,body->getId());
}

void LatticeExample::createActors(shared_ptr<MetaBody>& )
{
	shared_ptr<BoundingVolumeMetaEngine> boundingVolumeDispatcher	= shared_ptr<BoundingVolumeMetaEngine>(new BoundingVolumeMetaEngine);
	boundingVolumeDispatcher->add("MetaInteractingGeometry","AABB","MetaInteractingGeometry2AABB");

	shared_ptr<GeometricalModelMetaEngine> geometricalModelDispatcher	= shared_ptr<GeometricalModelMetaEngine>(new GeometricalModelMetaEngine);
	geometricalModelDispatcher->add("LatticeSetParameters","LatticeSetGeometry","LatticeSet2LatticeBeams");
	
//	strainRecorder = shared_ptr<StrainRecorder>(new StrainRecorder);
//	strainRecorder->outputFile 		= outputFile;
//	strainRecorder->interval 		= 10;
	
	rootBody->engines.clear();
	rootBody->engines.push_back(boundingVolumeDispatcher);
	rootBody->engines.push_back(shared_ptr<LatticeLaw>(new LatticeLaw));
	rootBody->engines.push_back(geometricalModelDispatcher);
//	rootBody->engines.push_back(strainRecorder);
	
	rootBody->initializers.clear();
	rootBody->initializers.push_back(boundingVolumeDispatcher);
	rootBody->initializers.push_back(geometricalModelDispatcher);
	
	if(useNonLocalModel)
	{
//		shared_ptr<NonLocalInitializer> nonLocalInitializer(new NonLocalInitializer);
//		nonLocalInitializer->range = nonLocalL_in_cellsizeUnit * cellsizeUnit_in_meters;
//		rootBody->initializers.push_back(nonLocalInitializer);
	}
}	

void LatticeExample::positionRootBody(shared_ptr<MetaBody>& rootBody)
{
	rootBody->isDynamic		= false;

	Quaternionr q;
	q.fromAxisAngle( Vector3r(0,0,1),0);
	shared_ptr<LatticeSetParameters> physics(new LatticeSetParameters);
	physics->se3			= Se3r(Vector3r(0,0,0),q);
	physics->beamGroupMask 		= beamGroupMask;
	physics->nodeGroupMask 		= nodeGroupMask;
	//physics->useBendTensileSoftening= useBendTensileSoftening;
	//physics->useStiffnessSoftening  = useStiffnessSoftening;
	
	shared_ptr<MetaInteractingGeometry> set(new MetaInteractingGeometry());
	
	set->diffuseColor		= Vector3f(0,0,1);

	shared_ptr<AABB> aabb(new AABB);
	aabb->diffuseColor		= Vector3r(0,0,1);

	shared_ptr<GeometricalModel> gm = dynamic_pointer_cast<GeometricalModel>(ClassFactory::instance().createShared("LatticeSetGeometry"));
	gm->diffuseColor 		= Vector3f(1,1,1);
	gm->wire 			= false;
	gm->visible 			= true;
	gm->shadowCaster 		= true;
	
	rootBody->interactionGeometry	= dynamic_pointer_cast<InteractingGeometry>(set);	
	rootBody->boundingVolume	= dynamic_pointer_cast<BoundingVolume>(aabb);
	rootBody->geometricalModel 	= gm;
	rootBody->physicalParameters 	= physics;
}
	
 
void LatticeExample::imposeTranslation(shared_ptr<MetaBody>& rootBody, Vector3r min, Vector3r max, Vector3r direction, Real displacement)
{
	shared_ptr<DisplacementEngine> translationCondition = shared_ptr<DisplacementEngine>(new DisplacementEngine);
 	translationCondition->displacement  = displacement;
	direction.normalize();
 	translationCondition->translationAxis = direction;
	
	rootBody->engines.push_back((rootBody->engines)[rootBody->engines.size()-1]);
	(rootBody->engines)[rootBody->engines.size()-2]=(rootBody->engines)[rootBody->engines.size()-3];
	(rootBody->engines)[rootBody->engines.size()-3]=(rootBody->engines)[rootBody->engines.size()-4];
	(rootBody->engines)[rootBody->engines.size()-4]=translationCondition;
	translationCondition->subscribedBodies.clear();
	
	BodyContainer::iterator bi    = rootBody->bodies->begin();
	BodyContainer::iterator biEnd = rootBody->bodies->end();
	for(  ; bi!=biEnd ; ++bi )
	{
		shared_ptr<Body> b = *bi;
	
		if( b->getGroupMask() & nodeGroupMask )
		{
			Vector3r pos = b->physicalParameters->se3.position;
			if(        pos[0] > min[0] 
				&& pos[1] > min[1] 
				&& pos[2] > min[2] 
				&& pos[0] < max[0] 
				&& pos[1] < max[1] 
				&& pos[2] < max[2] 
				&& (b->getGroupMask() & nodeGroupMask)
				)
			{
				b->isDynamic = false;
				b->geometricalModel->diffuseColor = Vector3f(0.3,0.3,0.3);
				translationCondition->subscribedBodies.push_back(b->getId());
			}
		}
	}
}

void LatticeExample::regionDelete(shared_ptr<MetaBody>& rootBody, Vector3r min, Vector3r max)
{
	vector<unsigned int> futureDeletes;
	
	BodyContainer::iterator bi    = rootBody->bodies->begin();
	BodyContainer::iterator biEnd = rootBody->bodies->end();
	for(  ; bi!=biEnd ; ++bi )
	{
		shared_ptr<Body> b = *bi;
	
		if( b->getGroupMask() & beamGroupMask )
		{
			Vector3r pos = b->physicalParameters->se3.position;
			if(        pos[0] > min[0] 
				&& pos[1] > min[1] 
				&& pos[2] > min[2] 
				&& pos[0] < max[0] 
				&& pos[1] < max[1] 
				&& pos[2] < max[2] 
				)
				futureDeletes.push_back( b->getId() );
		}
	}
	
	vector<unsigned int>::iterator vend = futureDeletes.end();
	for( vector<unsigned int>::iterator vsta = futureDeletes.begin() ; vsta != vend ; ++vsta)
		rootBody->bodies->erase(*vsta); 
}

void LatticeExample::nonDestroy(shared_ptr<MetaBody>& rootBody, Vector3r min, Vector3r max)
{
	vector<unsigned int> marked;
	
	BodyContainer::iterator bi    = rootBody->bodies->begin();
	BodyContainer::iterator biEnd = rootBody->bodies->end();
	for(  ; bi!=biEnd ; ++bi )
	{
		shared_ptr<Body> b = *bi;
	
		if( b->getGroupMask() & beamGroupMask )
		{
			Vector3r pos = b->physicalParameters->se3.position;
			if(        pos[0] > min[0] 
				&& pos[1] > min[1] 
				&& pos[2] > min[2] 
				&& pos[0] < max[0] 
				&& pos[1] < max[1] 
				&& pos[2] < max[2] 
				)
				marked.push_back( b->getId() );
		}
	}
	
	vector<unsigned int>::iterator vend = marked.end();
	for( vector<unsigned int>::iterator vsta = marked.begin() ; vsta != vend ; ++vsta)
	{
		LatticeBeamParameters* beam = static_cast<LatticeBeamParameters*>( ((*(rootBody->bodies))[*vsta])->physicalParameters.get());
		beam->criticalTensileStrain 	= 0.9;
		beam->criticalCompressiveStrain = 0.9;

		(*(rootBody->bodies))[beam->id1]->geometricalModel->diffuseColor = Vector3f(0.2,0.5,0.7);
		(*(rootBody->bodies))[beam->id2]->geometricalModel->diffuseColor = Vector3f(0.2,0.5,0.7);
	}
}




struct Circle
{
	float x,y,r;
};

bool overlaps(Circle& cc,std::vector<Circle>& c)
{
	std::vector<Circle>::iterator end=c.end();
	for(std::vector<Circle>::iterator i=c.begin();i!=end;++i)
	{
		float dist2 = std::pow(i->x - cc.x ,2)+std::pow(i->y - cc.y,2);
		float r2    = std::pow(i->r+cc.r ,2);
		if(dist2<r2)
			return true;
	}
	return false;
};

bool aggInside(Vector3r& a,Vector3r& b,std::vector<Circle>& c)
{ // checks if nodes 'a','b' are inside any of aggregates from list 'c'
	std::vector<Circle>::iterator end=c.end();
	for(std::vector<Circle>::iterator i=c.begin();i!=end;++i)
	{
		float dist2 = std::pow(i->x - a[0],2)+std::pow(i->y - a[1],2);
		float r2    = std::pow(i->r,2);
		if(dist2<r2)
		{
			dist2 = std::pow(i->x - b[0],2)+std::pow(i->y - b[1],2);
			if(dist2<r2)
				return true;
		}
	}
	return false;
}

float aggsAreas(std::vector<Circle>& c)
{
	float aggArea=0.0;
	std::vector<Circle>::iterator end=c.end();
	for(std::vector<Circle>::iterator i=c.begin();i!=end;++i)
		aggArea += 3.14159265358979323846*std::pow(i->r ,2);
	return aggArea;
}

// random
#include <boost/random/linear_congruential.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/normal_distribution.hpp>

void LatticeExample::addAggregates(shared_ptr<MetaBody>& rootBody)
{
	// first make a list of circles
	std::vector<Circle> c;
	float AGGREGATES_X=speciemen_size_in_meters[0];
	float AGGREGATES_Y=speciemen_size_in_meters[1];

	typedef boost::minstd_rand StdGenerator;
	static StdGenerator generator;
	static boost::variate_generator<StdGenerator&, boost::uniform_real<> >
		random1(generator, boost::uniform_real<>(0,1));
	static boost::variate_generator<StdGenerator&, boost::normal_distribution<> >
		randomN(generator, boost::normal_distribution<>(aggregateMeanRadius,aggregateSigmaRadius));

	std::cerr << "generating aggregates ... ";
	do
	{
		Circle cc;
		cc.x=random1()*AGGREGATES_X, cc.y=random1()*AGGREGATES_Y;
		cc.r=std::abs(randomN());
		for(int i=0 ; i<1000 ; ++i)
			if(overlaps(cc,c))
				cc.x=random1()*AGGREGATES_X, cc.y=random1()*AGGREGATES_Y;
			else
			{
				c.push_back(cc);
		//		std::cerr << cc.x << " " << cc.y << " " << cc.r << "\n";
				break;
			}
	}
	while(aggregatePercent/100.0 > aggsAreas(c)/(AGGREGATES_X*AGGREGATES_Y) );
	std::cerr << "done. " << c.size() << " aggregates generated, area: " << aggsAreas(c)/(AGGREGATES_X*AGGREGATES_Y) << "\n";

	{ // set different properties for beams that lie in an aggregate
		BodyContainer::iterator bi    = rootBody->bodies->begin();
		BodyContainer::iterator biEnd = rootBody->bodies->end();
		float all_bodies = rootBody->bodies->size();
		int current = 0;
		for(  ; bi!=biEnd ; ++bi )  // loop over all beams
		{
		//	if( ++current % 100 == 0 )
		//		cerr << "aggregate beams: " << current << " , " << ((static_cast<float>(current)/all_bodies)*100.0) << " %\n";
				
			Body* body = (*bi).get();
			if( ! ( body->getGroupMask() & beamGroupMask ) )
				continue; // skip non-beams

			LatticeBeamParameters* beam 	= static_cast<LatticeBeamParameters*>(body->physicalParameters.get());
			if(aggInside(		 (*(rootBody->bodies))[beam->id1]->physicalParameters->se3.position
						,(*(rootBody->bodies))[beam->id2]->physicalParameters->se3.position
						,c ))
			{
				beam->longitudalStiffness 	= agg_longStiffness_noUnit;
				beam->bendingStiffness 		= agg_bendStiffness_noUnit;
		
				(*(rootBody->bodies))[beam->id1]->geometricalModel->diffuseColor = Vector3f(0.0,0.0,0.0);
				(*(rootBody->bodies))[beam->id2]->geometricalModel->diffuseColor = Vector3f(0.0,0.0,0.0);
			}
		}
	}
}

/*************************************************************************
*  Copyright (C) 2004 by Janek Kozicki                                   *
*  cosurgi@xxxxxxxxxx                                                    *
*                                                                        *
*  This program is free software; it is licensed under the terms of the  *
*  GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/

#ifndef LATTICE_BOX_HPP
#define LATTICE_BOX_HPP 

#include <yade/yade-core/FileGenerator.hpp>
#include <yade/yade-lib-wm3-math/Vector3.hpp>
#include <yade/yade-package-common/BodyRedirectionVector.hpp>

class LatticeBeamParameters;
class StrainRecorder;

class LatticeExample : public FileGenerator
{
	private :
		int 		nodeGroupMask,beamGroupMask;
		
	// mesh generation	
		Vector3r 	 speciemen_size_in_meters 	// size
				,disorder_in_cellsizeUnit;	// s
		
		Real 		 maxLength_in_cellsizeUnit	// r
				,minAngle_betweenBeams_deg 	// a
				,cellsizeUnit_in_meters 	// g
				
				,crit_TensileStrain_percent	// E_min
				,crit_ComprStrain_percent	// E_max
				
				,longitudalStiffness_noUnit 	// k_l
				,bendingStiffness_noUnit	// k_b
				
				,nonLocalL_in_cellsizeUnit; 	// l
				
		bool 		 triangularBaseGrid
				,useNonLocalModel
				,useBendTensileSoftening
				,useStiffnessSoftening;
	
	// aggregates
		bool		 useAggregates;
		Real		 aggregatePercent
				,aggregateMeanRadius
				,aggregateSigmaRadius
				,agg_longStiffness_noUnit 	// k_l
				,agg_bendStiffness_noUnit;	// k_b
		void addAggregates(shared_ptr<MetaBody>& rootBody);
		
		
	// conditions
		Vector3r 	 region_A_min
			 	,region_A_max
			 	,direction_A
				 
			 	,region_B_min
			 	,region_B_max
			 	,direction_B
				 
			 	,region_C_min
			 	,region_C_max
			 	,direction_C
				 
			 	,region_D_min
			 	,region_D_max
			 	,direction_D;
				 
						
		Real		 displacement_A_meters
		 		,displacement_B_meters
		 		,displacement_C_meters
		 		,displacement_D_meters;
				
	// strain recorder
		std::vector< unsigned int > subscribedBodies;
		Real		 strainRecorder_xz_plane;
		shared_ptr<StrainRecorder> strainRecorder;
		std::string 	 outputFile;
		
	// delete beams regions
		Vector3r 	 regionDelete_A_min
			 	,regionDelete_A_max
			 	,regionDelete_B_min
			 	,regionDelete_B_max
				,regionDelete_C_min
			 	,regionDelete_C_max
				,regionDelete_D_min
			 	,regionDelete_D_max
				,regionDelete_E_min
			 	,regionDelete_E_max
				,regionDelete_F_min
			 	,regionDelete_F_max;

	// non destroy areas
		Vector3r	 nonDestroy_A_min
				,nonDestroy_A_max
				,nonDestroy_B_min
				,nonDestroy_B_max;
		
				 
		std::vector< std::vector< unsigned int > > connections; // which node is in touch with what beams.
				
	public : 
		LatticeExample();
		virtual ~LatticeExample();

		string generate();
	
		void createActors(shared_ptr<MetaBody>& rootBody);
		void positionRootBody(shared_ptr<MetaBody>& rootBody);
		bool createNode(shared_ptr<Body>& body, int i, int j, int k);
		void createBeam(shared_ptr<Body>& body, unsigned int i, unsigned int j);
		void calcBeamPositionOrientationLength(shared_ptr<Body>& body);
		void calcBeamAngles(Body* body, BodyContainer* bodies,InteractionContainer* ints);
		void calcAxisAngle(LatticeBeamParameters* beam, BodyContainer* bodies, unsigned int otherId,InteractionContainer* ints,unsigned int thisId);
		bool checkMinimumAngle(BodyRedirectionVector&,shared_ptr<Body>&);
		bool checkAngle( Vector3r , Vector3r& );
		void imposeTranslation(shared_ptr<MetaBody>& rootBody, Vector3r min, Vector3r max, Vector3r direction, Real velocity);
		void regionDelete(shared_ptr<MetaBody>& rootBody, Vector3r min, Vector3r max);
		void nonDestroy(shared_ptr<MetaBody>& rootBody, Vector3r min, Vector3r max);

		virtual void registerAttributes();
		REGISTER_CLASS_NAME(LatticeExample);
		REGISTER_BASE_CLASS_NAME(FileGenerator);

};

REGISTER_SERIALIZABLE(LatticeExample,false);

#endif // LATTICE_BOX_HPP