yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00783
[svn] r1573 - in trunk: extra/triangulation pkg/dem/Engine/DeusExMachina pkg/dem/Engine/StandAloneEngine
Author: chareyre
Date: 2008-11-17 16:41:43 +0100 (Mon, 17 Nov 2008)
New Revision: 1573
Modified:
trunk/extra/triangulation/Tesselation.cpp
trunk/extra/triangulation/TesselationWrapper.cpp
trunk/extra/triangulation/makefile
trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp
trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp
Log:
- Fix conflicts and prepare next implementation of isotropic compression
Modified: trunk/extra/triangulation/Tesselation.cpp
===================================================================
--- trunk/extra/triangulation/Tesselation.cpp 2008-11-17 10:54:38 UTC (rev 1572)
+++ trunk/extra/triangulation/Tesselation.cpp 2008-11-17 15:41:43 UTC (rev 1573)
@@ -1,13 +1,13 @@
#include "Tesselation.h"
#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
//Tesselation TESSELATION = Tesselation();
-
+
using namespace std;
-
+
Tesselation::Tesselation ( void )
{
-// std::cout << "Tesselation(void)" << std::endl;
+// std::cout << "Tesselation(void)" << std::endl;
Tri = new RTriangulation;
Tes = Tri;
computed = false;
@@ -22,7 +22,7 @@
Tesselation::Tesselation ( RTriangulation &T )
: Tri ( &T ), Tes ( &T ), computed ( true )
{
-// std::cout << "Tesselation(RTriangulation &T)" << std::endl;
+// std::cout << "Tesselation(RTriangulation &T)" << std::endl;
Compute();
}
@@ -42,20 +42,20 @@
// bool Tesselation::insert ( Real x, Real y, Real z, Real rad, unsigned int id, bool isFictious )
// {
-// Vertex_handle Vh;
-// Vh = Tri->insert ( Sphere ( Point ( x,y,z ),pow ( rad,2 ) ) );
-// if ( Vh!=NULL )
-// {
-// Vh->info() = id;
-// Vh->info().isFictious = isFictious;
-// if ( !isFictious ) max_id = std::max ( max_id, id );
-// }
-// else return false;
-// return true;
+// Vertex_handle Vh;
+// Vh = Tri->insert ( Sphere ( Point ( x,y,z ),pow ( rad,2 ) ) );
+// if ( Vh!=NULL )
+// {
+// Vh->info() = id;
+// Vh->info().isFictious = isFictious;
+// if ( !isFictious ) max_id = std::max ( max_id, id );
+// }
+// else return false;
+// return true;
// }
-Vertex_handle Tesselation::insert(Real x, Real y, Real z, Real rad, unsigned int id, bool isFictious)
+Vertex_handle Tesselation::insert ( Real x, Real y, Real z, Real rad, unsigned int id, bool isFictious )
{
Vertex_handle Vh;
@@ -66,12 +66,12 @@
Vh->info().isFictious = isFictious;
if ( !isFictious ) max_id = std::max ( max_id, id );
}
- else cout << id << " : Vh==NULL!!" << endl;
+ else cout << id << " : Vh==NULL!!" << endl;
return Vh;
}
-Vertex_handle Tesselation::move(Real x, Real y, Real z, Real rad, unsigned int id)
+Vertex_handle Tesselation::move ( Real x, Real y, Real z, Real rad, unsigned int id )
{
bool fictious = vertexHandles[id]->info().isFictious;
Vertex_handle Vh;
@@ -80,13 +80,14 @@
{
vertexHandles[id] = Vh;
Vh->info() = id;
- Vh->info().isFictious = fictious;
- } else cerr << "Vh==NULL" << " id=" << id << " Point=" << Point ( x,y,z ) << " rad=" << rad << endl;
+ Vh->info().isFictious = fictious;
+ }
+ else cerr << "Vh==NULL" << " id=" << id << " Point=" << Point ( x,y,z ) << " rad=" << rad << endl;
return Vh;
}
-
+
bool Tesselation::redirect ( void )
{
if ( !redirected )
@@ -159,11 +160,8 @@
void Tesselation::Compute ()
{
-<<<<<<< .mine
-// std::cout << "Tesselation::Compute ()" << std::endl;
-=======
- std::cerr << "Tesselation::Compute ()" << std::endl;
->>>>>>> .r1529
+
+// std::cout << "Tesselation::Compute ()" << std::endl;
Finite_cells_iterator cell_end = Tri->finite_cells_end();
for ( Finite_cells_iterator cell = Tri->finite_cells_begin(); cell != cell_end; cell++ )
@@ -185,11 +183,11 @@
S3.point().x(), S3.point().y(), S3.point().z(), S3.weight(),
x, y, z );
cell->info() =Point ( x,y,z );
- //cout << "voronoi cell : " << cell->vertex(0)->info().id() << " "
- // << cell->vertex(1)->info().id() << " "
- // << cell->vertex(2)->info().id() << " "
- // << cell->vertex(3)->info().id() << "(center : " << (Point) cell->info() << ")" << endl;
-}
+ //cout << "voronoi cell : " << cell->vertex(0)->info().id() << " "
+ // << cell->vertex(1)->info().id() << " "
+ // << cell->vertex(2)->info().id() << " "
+ // << cell->vertex(3)->info().id() << "(center : " << (Point) cell->info() << ")" << endl;
+ }
computed = true;
@@ -400,30 +398,30 @@
}
cell0=cell2++;
Cell_circulator cell1=cell2++;
- //std::cout << "edge : " << ed_it->first->vertex ( ed_it->second )->info().id() << "-" << ed_it->first->vertex ( ed_it->third )->info().id() << std::endl;
+ //std::cout << "edge : " << ed_it->first->vertex ( ed_it->second )->info().id() << "-" << ed_it->first->vertex ( ed_it->third )->info().id() << std::endl;
bool isFictious1 = ( ed_it->first )->vertex ( ed_it->second )->info().isFictious;
bool isFictious2 = ( ed_it->first )->vertex ( ed_it->third )->info().isFictious;
Real r;
-
- //cout << "cell0 : " << cell0->vertex(0)->info().id() << " "
- // << cell0->vertex(1)->info().id() << " "
- // << cell0->vertex(2)->info().id() << " "
- // << cell0->vertex(3)->info().id() << "(center : " << (Point) cell0->info() << ")" << endl;
+ //cout << "cell0 : " << cell0->vertex(0)->info().id() << " "
+ // << cell0->vertex(1)->info().id() << " "
+ // << cell0->vertex(2)->info().id() << " "
+ // << cell0->vertex(3)->info().id() << "(center : " << (Point) cell0->info() << ")" << endl;
+
while ( cell2!=cell0 )
- {
+ {
if ( !Tri->is_infinite ( cell1 ) && !Tri->is_infinite ( cell2 ) )
{
- // cout << "cell1 : " << cell1->vertex(0)->info().id() << " "
- // << cell1->vertex(1)->info().id() << " "
- // << cell1->vertex(2)->info().id() << " "
- // << cell1->vertex(3)->info().id() << "(center : " << (Point) cell1->info() << ")" << endl;
- // cout << "cell2 : " << cell2->vertex(0)->info().id() << " "
- // << cell2->vertex(1)->info().id() << " "
- // << cell2->vertex(2)->info().id() << " "
- // << cell2->vertex(3)->info().id() << "(center : " << (Point) cell2->info() << ")" << endl;
-
- //std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->second )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
+ // cout << "cell1 : " << cell1->vertex(0)->info().id() << " "
+ // << cell1->vertex(1)->info().id() << " "
+ // << cell1->vertex(2)->info().id() << " "
+ // << cell1->vertex(3)->info().id() << "(center : " << (Point) cell1->info() << ")" << endl;
+ // cout << "cell2 : " << cell2->vertex(0)->info().id() << " "
+ // << cell2->vertex(1)->info().id() << " "
+ // << cell2->vertex(2)->info().id() << " "
+ // << cell2->vertex(3)->info().id() << "(center : " << (Point) cell2->info() << ")" << endl;
+
+ //std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->second )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
if ( !isFictious1 )
{
r = std::abs ( ( Tetraedre ( ed_it->first->vertex ( ed_it->second )->point(), cell0->info(), cell1->info(), cell2->info() ) ).volume() );
@@ -431,7 +429,7 @@
( ed_it->first )->vertex ( ed_it->second )->info().v() += r;
TotalFiniteVoronoiVolume+=r;
}
- //std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->third )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
+ //std::cout << "assign tetra : (" << ed_it->first->vertex ( ed_it->third )->point() << "),(" << cell0->info() << "),(" << cell1->info() << "),(" <<cell2->info() << ")" << std::endl;
if ( !isFictious2 )
{
r = std::abs ( ( Tetraedre ( ed_it->first->vertex ( ed_it->third )->point(), cell0->info(), cell1->info(), cell2->info() ) ).volume() );
@@ -474,11 +472,12 @@
redirect();
-// cout << "TotalVolume : " << TotalFiniteVoronoiVolume << endl;
+// cout << "TotalVolume : " << TotalFiniteVoronoiVolume << endl;
}
void Tesselation::ComputePorosity ( void ) //WARNING : This function will erase real volumes of cells
-{ // and replace it with porosity
+{
+ // and replace it with porosity
ComputeVolumes();
//Real rr=0;
Finite_vertices_iterator vertices_end = Tri->finite_vertices_end ();
Modified: trunk/extra/triangulation/TesselationWrapper.cpp
===================================================================
--- trunk/extra/triangulation/TesselationWrapper.cpp 2008-11-17 10:54:38 UTC (rev 1572)
+++ trunk/extra/triangulation/TesselationWrapper.cpp 2008-11-17 15:41:43 UTC (rev 1573)
@@ -57,11 +57,9 @@
facet_it = Tes->Triangulation().finite_edges_end ();
}
-<<<<<<< .mine
+
-=======
-
void TesselationWrapper::clear2(void) //for testing purpose
{
Tes->Clear();
@@ -76,7 +74,7 @@
}
->>>>>>> .r1529
+
double TesselationWrapper::Volume( unsigned int id )
{
return Tes->Volume(id);
Modified: trunk/extra/triangulation/makefile
===================================================================
--- trunk/extra/triangulation/makefile 2008-11-17 10:54:38 UTC (rev 1572)
+++ trunk/extra/triangulation/makefile 2008-11-17 15:41:43 UTC (rev 1573)
@@ -83,18 +83,12 @@
libTesselationWrapper.a : TesselationWrapper$(OBJ_EXT) \
Tenseur3$(OBJ_EXT) \
Tesselation$(OBJ_EXT) \
-<<<<<<< .mine
- RegularTriangulation$(OBJ_EXT)
- ar cr libTesselationWrapper.a \
- TesselationWrapper$(OBJ_EXT) \
-=======
RegularTriangulation$(OBJ_EXT) \
TriaxialState$(OBJ_EXT) \
KinematicLocalisationAnalyser$(OBJ_EXT)
ar cr libTesselationWrapper.a \
TesselationWrapper$(OBJ_EXT) \
Tenseur3$(OBJ_EXT) \
->>>>>>> .r1529
Tesselation$(OBJ_EXT) \
RegularTriangulation$(OBJ_EXT) \
TriaxialState$(OBJ_EXT) \
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp 2008-11-17 10:54:38 UTC (rev 1572)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialCompressionEngine.cpp 2008-11-17 15:41:43 UTC (rev 1573)
@@ -51,14 +51,14 @@
previousSigmaIso=sigma_iso;
frictionAngleDegree = -1;
epsilonMax = 0.5;
+<<<<<<< .mine
isotropicCompaction=false;
spheresVolume=0;
boxVolume=0;
- calculatedPorosity=1.1;
+// calculatedPorosity=1.1;
-
}
TriaxialCompressionEngine::~TriaxialCompressionEngine()
@@ -184,11 +184,13 @@
doStateTransition (ncb, STATE_LIMBO );
}
}
- else if ( calculatedPorosity<=fixedPorosity && currentState==STATE_FIXED_POROSITY_COMPACTION )
+
+ else if ( porosity<=fixedPorosity && currentState==STATE_FIXED_POROSITY_COMPACTION )
{
- Omega::instance().stopSimulationLoop();
- return;
+ Omega::instance().stopSimulationLoop();
+ return;
}
+
#if 0
//This is a hack in order to allow subsequent run without activating compression - like for the YADE-COMSOL coupling
if ( !compressionActivated )
@@ -222,21 +224,22 @@
LOG_INFO ( "First run, will initialize!" );
/* FIXME: are these three if's mutually exclusive and are partition of all possibilities? */
//sigma_iso was changed, we need to rerun compaction
+
if ( (sigmaIsoCompaction!=previousSigmaIso || currentState==STATE_UNINITIALIZED || currentState== STATE_LIMBO) && currentState!=STATE_TRIAX_LOADING && isotropicCompaction == false) doStateTransition (ncb, STATE_ISO_COMPACTION );
if ( previousState==STATE_LIMBO && currentState==STATE_TRIAX_LOADING && isotropicCompaction == false ) doStateTransition (ncb, STATE_TRIAX_LOADING );
- if ( fixedPorosity<1 && currentState==STATE_UNINITIALIZED && isotropicCompaction!=false )
- {
- doStateTransition (ncb, STATE_FIXED_POROSITY_COMPACTION );
- FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
- if(!b->isDynamic) continue;
- const shared_ptr<Sphere>& sphere=YADE_PTR_CAST<Sphere> ( b->geometricalModel );
- spheresVolume += (4./3.)*Mathr::PI*pow( sphere->radius, 3 );
- }
- previousState=currentState;
- previousSigmaIso=sigma_iso;
- firstRun=false; // change this only _after_ state transitions
- }
+ //if ( fixedPorosity<1 && currentState==STATE_UNINITIALIZED && isotropicCompaction!=false )
+ //{
+ // doStateTransition (ncb, STATE_FIXED_POROSITY_COMPACTION );
+// FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
+// if(!b->isDynamic) continue;
+// const shared_ptr<Sphere>& sphere=YADE_PTR_CAST<Sphere> ( b->geometricalModel );
+// spheresVolume += (4./3.)*Mathr::PI*pow( sphere->radius, 3 );
+// }
+ previousState=currentState;
+ previousSigmaIso=sigma_iso;
+ firstRun=false; // change this only _after_ state transitions
}
+
if ( saveSimulation )
{
string fileName = "./"+ Key + "_" + Phase1End + "_" +
@@ -288,6 +291,7 @@
Omega::instance().stopSimulationLoop();
}
}
+
if ( currentState==STATE_FIXED_POROSITY_COMPACTION )
{
if ( Omega::instance().getCurrentIteration() % 100 == 0 )
@@ -311,13 +315,6 @@
p->se3.position += 0.5*translationSpeed*width*translationAxisx*dt;
p = static_cast<PhysicalParameters*> ( Body::byId ( wall_right_id )->physicalParameters.get() );
p->se3.position -= 0.5*translationSpeed*width*translationAxisx*dt;
-
-
- boxVolume=height*width*depth;
-
- calculatedPorosity=1-spheresVolume/boxVolume;
-
-
}
}
Modified: trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2008-11-17 10:54:38 UTC (rev 1572)
+++ trunk/pkg/dem/Engine/DeusExMachina/TriaxialStressController.cpp 2008-11-17 15:41:43 UTC (rev 1573)
@@ -27,15 +27,17 @@
//StiffnessMatrixClassIndex = actionParameterStiffnessMatrix->getClassIndex();
shared_ptr<Force> tmpF(new Force);
ForceClassIndex=tmpF->getClassIndex();
-
+ firstRun = true;
+
previousStress = 0;
previousMultiplier = 1;
-
+
stiffnessUpdateInterval =10;
radiusControlInterval =10;
computeStressStrainInterval = 10;
wallDamping = 0.25;
//force = Vector3r::ZERO;
+ stiffness.resize(6);
for (int i=0; i<6; ++i)
{
wall_id[i] = 0;
@@ -68,6 +70,9 @@
height = 0;
width = 0;
depth = 0;
+ spheresVolume=0;
+ boxVolume=0;
+ porosity=0;
height0 = 0;
width0 = 0;
depth0 = 0;
@@ -93,7 +98,7 @@
//REGISTER_ATTRIBUTE(UnbalancedForce);
- //REGISTER_ATTRIBUTE(stiffness);
+ REGISTER_ATTRIBUTE(stiffness);
REGISTER_ATTRIBUTE(wall_bottom_id);
REGISTER_ATTRIBUTE(wall_top_id);
REGISTER_ATTRIBUTE(wall_left_id);
@@ -193,8 +198,11 @@
void TriaxialStressController::applyCondition(MetaBody* ncb)
{
//cerr << "TriaxialStressController::applyCondition" << endl;
-
+
+
shared_ptr<BodyContainer>& bodies = ncb->bodies;
+
+
if(thickness<=0) thickness=YADE_PTR_CAST<InteractingBox>(Body::byId(wall_bottom_id,ncb)->interactingGeometry)->extents.Y();
@@ -208,8 +216,28 @@
height = p_top->se3.position.Y() - p_bottom->se3.position.Y() - thickness;
width = p_right->se3.position.X() - p_left->se3.position.X() - thickness;
depth = p_front->se3.position.Z() - p_back->se3.position.Z() - thickness;
+
+ boxVolume = height * width * depth;
+
+ if (firstRun) {
+ BodyContainer::iterator bi = ncb->bodies->begin();
+ BodyContainer::iterator biEnd = ncb->bodies->end();
+ spheresVolume = 0;
+ for ( ; bi!=biEnd; ++bi )
+ {
+ const shared_ptr<Body>& b = *bi;
+ if ( b->isDynamic )
+ {
+ const shared_ptr<Sphere>& sphere =
+ YADE_PTR_CAST<Sphere> ( b->geometricalModel );
+ spheresVolume += 1.3333333*Mathr::PI*pow ( sphere->radius, 3 );
+ }
+ }
+
+ firstRun = false;
+ }
+ //porosity = ( boxVolume - spheresVolume ) /boxVolume;
-
position_top = p_top->se3.position.Y();
position_bottom = p_bottom->se3.position.Y();
position_right = p_right->se3.position.X();
@@ -221,7 +249,7 @@
// must be done _after_ height, width, depth have been calculated
//Update stiffness only if it has been computed by StiffnessCounter (see "stiffnessUpdateInterval")
- if (Omega::instance().getCurrentIteration() % stiffnessUpdateInterval == 0 || Omega::instance().getCurrentIteration()<1000) updateStiffness(ncb);
+ if (Omega::instance().getCurrentIteration() % stiffnessUpdateInterval == 0 || Omega::instance().getCurrentIteration()<100) updateStiffness(ncb);
bool isARadiusControlIteration = (Omega::instance().getCurrentIteration() % radiusControlInterval == 0);
@@ -233,9 +261,11 @@
Vector3r wallForce (0, sigma_iso*width*depth, 0);
if (wall_bottom_activated) controlExternalStress(wall_bottom, ncb, -wallForce, p_bottom, max_vel);
if (wall_top_activated) controlExternalStress(wall_top, ncb, wallForce, p_top, max_vel);
+
wallForce = Vector3r(sigma_iso*height*depth, 0, 0);
if (wall_left_activated) controlExternalStress(wall_left, ncb, -wallForce, p_left, max_vel*width/height);
if (wall_right_activated) controlExternalStress(wall_right, ncb, wallForce, p_right, max_vel*width/height);
+
wallForce = Vector3r(0, 0, sigma_iso*height*width);
if (wall_back_activated) controlExternalStress(wall_back, ncb, -wallForce, p_back, max_vel*depth/height);
if (wall_front_activated) controlExternalStress(wall_front, ncb, wallForce, p_front, max_vel*depth/height);
@@ -316,39 +346,40 @@
meanStress/=6.;
}
-void TriaxialStressController::controlInternalStress(MetaBody* ncb, Real multiplier)
+void TriaxialStressController::controlInternalStress ( MetaBody* ncb, Real multiplier )
{
- BodyContainer::iterator bi = ncb->bodies->begin();
- BodyContainer::iterator biEnd = ncb->bodies->end();
- //cerr << "meanstress = "radius = " << endl;
- //cerr << "bouclesurBodies" << endl;
- for ( ; bi!=biEnd ; ++bi )
- {
- if ((*bi)->isDynamic)
+ spheresVolume *= pow ( multiplier,3 );
+ BodyContainer::iterator bi = ncb->bodies->begin();
+ BodyContainer::iterator biEnd = ncb->bodies->end();
+ //cerr << "meanstress = "radius = " << endl;
+ //cerr << "bouclesurBodies" << endl;
+ for ( ; bi!=biEnd ; ++bi )
+ {
+ if ( ( *bi )->isDynamic )
{
- (static_cast<InteractingSphere*> ((*bi)->interactingGeometry.get()))->radius *= multiplier;
- (static_cast<Sphere*>((*bi)->geometricalModel.get()))->radius *= multiplier;
- (static_cast<ParticleParameters*>((*bi)->physicalParameters.get()))->mass *= pow(multiplier,3);
- (static_cast<RigidBodyParameters*>((*bi)->physicalParameters.get()))->inertia *= pow(multiplier,5);
-
+ ( static_cast<InteractingSphere*> ( ( *bi )->interactingGeometry.get() ) )->radius *= multiplier;
+ ( static_cast<Sphere*> ( ( *bi )->geometricalModel.get() ) )->radius *= multiplier;
+ ( static_cast<ParticleParameters*> ( ( *bi )->physicalParameters.get() ) )->mass *= pow ( multiplier,3 );
+ ( static_cast<RigidBodyParameters*> ( ( *bi )->physicalParameters.get() ) )->inertia *= pow ( multiplier,5 );
+
}
}
// << "bouclesurInteraction" << endl;
InteractionContainer::iterator ii = ncb->transientInteractions->begin();
InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
- for ( ; ii!=iiEnd ; ++ii )
+ for ( ; ii!=iiEnd ; ++ii )
{
- if ((*ii)->isReal)
+ if ( ( *ii )->isReal )
{
- SpheresContactGeometry* contact = static_cast<SpheresContactGeometry*> ((*ii)->interactionGeometry.get());
- // if ((*(ncb->bodies))[(*ii)->getId1()]->isDynamic)
- // contact->radius1 *= multiplier;
- // if ((*(ncb->bodies))[(*ii)->getId2()]->isDynamic)
- // contact->radius2 *= multiplier;
- if ((*(ncb->bodies))[(*ii)->getId1()]->isDynamic)
- contact->radius1 = static_cast<InteractingSphere*> ((*(ncb->bodies))[(*ii)->getId1()]->interactingGeometry.get())->radius;
- if ((*(ncb->bodies))[(*ii)->getId2()]->isDynamic)
- contact->radius2 = static_cast<InteractingSphere*> ((*(ncb->bodies))[(*ii)->getId2()]->interactingGeometry.get())->radius;
+ SpheresContactGeometry* contact = static_cast<SpheresContactGeometry*> ( ( *ii )->interactionGeometry.get() );
+ // if ((*(ncb->bodies))[(*ii)->getId1()]->isDynamic)
+ // contact->radius1 *= multiplier;
+ // if ((*(ncb->bodies))[(*ii)->getId2()]->isDynamic)
+ // contact->radius2 *= multiplier;
+ if ( ( * ( ncb->bodies ) ) [ ( *ii )->getId1() ]->isDynamic )
+ contact->radius1 = static_cast<InteractingSphere*> ( ( * ( ncb->bodies ) ) [ ( *ii )->getId1() ]->interactingGeometry.get() )->radius;
+ if ( ( * ( ncb->bodies ) ) [ ( *ii )->getId2() ]->isDynamic )
+ contact->radius2 = static_cast<InteractingSphere*> ( ( * ( ncb->bodies ) ) [ ( *ii )->getId2() ]->interactingGeometry.get() )->radius;
}
}
}
Modified: trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp 2008-11-17 10:54:38 UTC (rev 1572)
+++ trunk/pkg/dem/Engine/StandAloneEngine/VolumicContactLaw.cpp 2008-11-17 15:41:43 UTC (rev 1573)
@@ -19,13 +19,10 @@
#include <yade/pkg-common/InteractingSphere.hpp>
-
#include "VolumicContactLaw.hpp"
-
#include<yade/extra/TesselationWrapper.h>
+#include <time.h>
-
-
VolumicContactLaw::VolumicContactLaw() : InteractionSolver() , actionForce(new Force) , actionMomentum(new Momentum)
{
sdecGroupMask=1;
@@ -44,12 +41,79 @@
}
+
+void VolumicContactLaw::speedTest(MetaBody* ncb)
+{
+//BEGIN SPEED TEST
+ shared_ptr<BodyContainer>& bodies = ncb->bodies;
+ TesselationWrapper T;
+ BodyContainer::iterator biBegin = bodies->begin();
+ BodyContainer::iterator biEnd = bodies->end();
+ BodyContainer::iterator bi = biBegin;
+ for( ; bi!=biEnd ; ++bi )
+ {
+ if ((*bi)->isDynamic) {//means "is it a sphere (not a wall)"
+ const InteractingSphere* s = YADE_CAST<InteractingSphere*>((*bi)->interactingGeometry.get());
+ const RigidBodyParameters* p = YADE_CAST<RigidBodyParameters*>((*bi)->physicalParameters.get());
+ T.checkMinMax(p->se3.position[0],p->se3.position[1],p->se3.position[2], s->radius);
+ }
+ }
+
+ clock_t T1 = clock();
+
+ for (int j=0; j<30; j++)
+ {
+ T.clear2();
+ T.bounded = false;
+ for( bi = biBegin; bi!=biEnd ; ++bi )
+ {
+ if ((*bi)->isDynamic) {//means "is it a sphere (not a wall)"
+ const InteractingSphere* s = YADE_CAST<InteractingSphere*>((*bi)->interactingGeometry.get());
+ const RigidBodyParameters* p = YADE_CAST<RigidBodyParameters*>((*bi)->physicalParameters.get());
+ T.insert(p->se3.position[0],p->se3.position[1],p->se3.position[2], s->radius, (*bi)->getId());
+ }
+ }
+ //T.AddBoundingPlanes();
+ }
+
+
+ cerr << "Bouding planes apres : time = "<< difftime(clock(), T1)/ CLOCKS_PER_SEC << endl;
+ T1 = clock();
+
+
+ for (int j=0; j<30; j++)
+ {
+ T.clear2();
+ T.bounded = false;
+ //T.AddBoundingPlanes();
+ for( bi = biBegin; bi!=biEnd ; ++bi )
+ {
+ if ((*bi)->isDynamic) {//means "is it a sphere (not a wall)"
+ const InteractingSphere* s = YADE_CAST<InteractingSphere*>((*bi)->interactingGeometry.get());
+ const RigidBodyParameters* p = YADE_CAST<RigidBodyParameters*>((*bi)->physicalParameters.get());
+ T.insert(p->se3.position[0],p->se3.position[1],p->se3.position[2], s->radius, (*bi)->getId());
+ }
+
+ }
+
+ }
+
+ cerr << "Bouding planes avant : time = "<< difftime(clock(), T1)/ CLOCKS_PER_SEC << endl;
+//END SPEED TEST
+
+
+
+
+}
+
+
void VolumicContactLaw::action(MetaBody* ncb)
{
shared_ptr<BodyContainer>& bodies = ncb->bodies;
Real dt = Omega::instance().getTimeStep();
+ speedTest(ncb);
//BEGIN VORONOI TESSELATION
TesselationWrapper T1;
@@ -139,7 +203,6 @@
}
}
//ENDOF VORONOI TESSELATION
-
/// Non Permanents Links ///
