yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01333
[svn] r1813 - in trunk/pkg/dem: . Engine/DeusExMachina Engine/StandAloneEngine PreProcessor
Author: chareyre
Date: 2009-06-24 19:04:39 +0200 (Wed, 24 Jun 2009)
New Revision: 1813
Added:
trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.cpp
trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.hpp
Modified:
trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp
trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp
trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
trunk/pkg/dem/PreProcessor/TriaxialTestWater.cpp
trunk/pkg/dem/SConscript
Log:
1. Update capillary files and add a new engine for them.
2. Prepare for removing old useless classes (typically some variants of engines with "water" at the end)
Modified: trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp 2009-06-24 16:28:58 UTC (rev 1812)
+++ trunk/pkg/dem/Engine/DeusExMachina/CapillaryStressRecorder.cpp 2009-06-24 17:04:39 UTC (rev 1813)
@@ -43,7 +43,10 @@
{
if(deserializing)
{
- ofile.open(outputFile.c_str());
+ //bool file_exists = std::ifstream (outputFile.c_str()); //if file does not exist, we will write colums titles
+ ofile.open(outputFile.c_str(), std::ios::app);
+ //if (!file_exists) ofile<<"iteration s11 s22 s33 e11 e22 e33 unb_force porosity kineticE" << endl;
+
}
}
Modified: trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp 2009-06-24 16:28:58 UTC (rev 1812)
+++ trunk/pkg/dem/Engine/DeusExMachina/ContactStressRecorder.cpp 2009-06-24 17:04:39 UTC (rev 1813)
@@ -9,6 +9,7 @@
#include "ContactStressRecorder.hpp"
#include <yade/pkg-common/RigidBodyParameters.hpp>
#include <yade/pkg-common/ParticleParameters.hpp>
+//#include <yade/pkg-common/Force.hpp>
#include <yade/pkg-common/Sphere.hpp>
#include <yade/pkg-dem/BodyMacroParameters.hpp>
#include <yade/pkg-dem/ElasticContactLaw.hpp>
@@ -22,35 +23,34 @@
#include <yade/core/MetaBody.hpp>
#include <boost/lexical_cast.hpp>
-CREATE_LOGGER(ContactStressRecorder);
+CREATE_LOGGER ( ContactStressRecorder );
-ContactStressRecorder::ContactStressRecorder () : DataRecorder()
-
+ContactStressRecorder::ContactStressRecorder () : DataRecorder()/*, actionForce ( new Force )*/
{
outputFile = "";
interval = 1;
-
+
height = 0;
width = 0;
depth = 0;
thickness = 0;
- upperCorner = Vector3r(0,0,0);
- lowerCorner = Vector3r(0,0,0);
-
- sphere_ptr = shared_ptr<GeometricalModel> (new Sphere);
+ upperCorner = Vector3r ( 0,0,0 );
+ lowerCorner = Vector3r ( 0,0,0 );
+
+ sphere_ptr = shared_ptr<GeometricalModel> ( new Sphere );
SpheresClassIndex = sphere_ptr->getClassIndex();
-
+
//triaxCompEng = NULL;
//sampleCapPressEng = new SampleCapillaryPressureEngine;
-
+
}
-void ContactStressRecorder::postProcessAttributes(bool deserializing)
+void ContactStressRecorder::postProcessAttributes ( bool deserializing )
{
- if(deserializing)
+ if ( deserializing )
{
- ofile.open(outputFile.c_str());
+ ofile.open ( outputFile.c_str() );
}
}
@@ -58,36 +58,36 @@
void ContactStressRecorder::registerAttributes()
{
DataRecorder::registerAttributes();
- REGISTER_ATTRIBUTE(outputFile);
- REGISTER_ATTRIBUTE(interval);
-
- REGISTER_ATTRIBUTE(wall_bottom_id);
- REGISTER_ATTRIBUTE(wall_top_id);
- REGISTER_ATTRIBUTE(wall_left_id);
- REGISTER_ATTRIBUTE(wall_right_id);
- REGISTER_ATTRIBUTE(wall_front_id);
- REGISTER_ATTRIBUTE(wall_back_id);
-
- REGISTER_ATTRIBUTE(height);
- REGISTER_ATTRIBUTE(width);
- REGISTER_ATTRIBUTE(depth);
- REGISTER_ATTRIBUTE(thickness);
- REGISTER_ATTRIBUTE(upperCorner);
- REGISTER_ATTRIBUTE(lowerCorner);
+ REGISTER_ATTRIBUTE ( outputFile );
+ REGISTER_ATTRIBUTE ( interval );
+ REGISTER_ATTRIBUTE ( wall_bottom_id );
+ REGISTER_ATTRIBUTE ( wall_top_id );
+ REGISTER_ATTRIBUTE ( wall_left_id );
+ REGISTER_ATTRIBUTE ( wall_right_id );
+ REGISTER_ATTRIBUTE ( wall_front_id );
+ REGISTER_ATTRIBUTE ( wall_back_id );
+
+ REGISTER_ATTRIBUTE ( height );
+ REGISTER_ATTRIBUTE ( width );
+ REGISTER_ATTRIBUTE ( depth );
+ REGISTER_ATTRIBUTE ( thickness );
+ REGISTER_ATTRIBUTE ( upperCorner );
+ REGISTER_ATTRIBUTE ( lowerCorner );
+
}
bool ContactStressRecorder::isActivated()
{
- return ((Omega::instance().getCurrentIteration() % interval == 0) && (ofile));
+ return ( ( Omega::instance().getCurrentIteration() % interval == 0 ) && ( ofile ) );
}
-void ContactStressRecorder::action(MetaBody * ncb)
+void ContactStressRecorder::action ( MetaBody * ncb )
{
shared_ptr<BodyContainer>& bodies = ncb->bodies;
-
+
if ( !triaxCompEng )
{
vector<shared_ptr<Engine> >::iterator itFirst = ncb->engines.begin();
@@ -103,259 +103,253 @@
}
if ( !triaxCompEng ) LOG_DEBUG ( "stress controller engine NOT found" );
}
-
+
Real f1_el_x=0, f1_el_y=0, f1_el_z=0, x1=0, y1=0, z1=0, x2=0, y2=0, z2=0;
-
- Real sig11_el=0, sig22_el=0, sig33_el=0, sig12_el=0, sig13_el=0,
- sig23_el=0;
+
+ Real sig11_el=0, sig22_el=0, sig33_el=0, sig12_el=0, sig13_el=0, sig23_el=0;
//, Vwater = 0,
Real kinematicE = 0;
-
+
InteractionContainer::iterator ii = ncb->transientInteractions->begin();
- InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
-
- Real j = 0;
- Real FT[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
-
- for( ; ii!=iiEnd ; ++ii )
- {
- if ((*ii)->isReal())
- {
- const shared_ptr<Interaction>& interaction = *ii;
-
- unsigned int id1 = interaction -> getId1();
+ InteractionContainer::iterator iiEnd = ncb->transientInteractions->end();
+
+ Real j = 0;
+ Real FT[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
+
+ for ( ; ii!=iiEnd ; ++ii )
+ {
+ if ( ( *ii )->isReal() )
+ {
+ const shared_ptr<Interaction>& interaction = *ii;
+
+ unsigned int id1 = interaction -> getId1();
unsigned int id2 = interaction -> getId2();
-
- SpheresContactGeometry* currentContactGeometry =
- static_cast<SpheresContactGeometry*>(interaction->interactionGeometry.
- get());
-
+
+ SpheresContactGeometry* currentContactGeometry =
+ static_cast<SpheresContactGeometry*> ( interaction->interactionGeometry.get() );
+
ElasticContactInteraction* currentContactPhysics =
- static_cast<ElasticContactInteraction*>
- (interaction->interactionPhysics.get());
-
+ static_cast<ElasticContactInteraction*> ( interaction->interactionPhysics.get() );
+
Real fn = currentContactPhysics->normalForce.Length();
- if (fn!=0)
-
- //if (currentContactGeometry->penetrationDepth >= 0)
-
- {
- j=j+1;
-
- Vector3r fel =
-currentContactPhysics->normalForce + currentContactPhysics->shearForce;
-
- f1_el_x=fel[0];
- f1_el_y=fel[1];
- f1_el_z=fel[2];
-
- int geometryIndex1 =
- (*bodies)[id1]->geometricalModel->getClassIndex();
- int geometryIndex2 =
- (*bodies)[id2]->geometricalModel->getClassIndex();
-
- if (geometryIndex1 == geometryIndex2)
-
- { BodyMacroParameters* de1 =
-static_cast<BodyMacroParameters*>((*bodies)[id1]->physicalParameters.get());
- x1 = de1->se3.position[0];
- y1 = de1->se3.position[1];
- z1 = de1->se3.position[2];
+ if ( fn!=0 )
+ //if (currentContactGeometry->penetrationDepth >= 0)
- BodyMacroParameters* de2 =
-static_cast<BodyMacroParameters*>((*bodies)[id2]->physicalParameters.get());
- x2 = de2->se3.position[0];
- y2 = de2->se3.position[1];
- z2 = de2->se3.position[2];
-
- ///Calcul des contraintes elastiques spheres/spheres
-
- sig11_el = sig11_el + f1_el_x*(x2 - x1);
- sig22_el = sig22_el + f1_el_y*(y2 - y1);
- sig33_el = sig33_el + f1_el_z*(z2 - z1);
- sig12_el = sig12_el + f1_el_x*(y2 - y1);
- sig13_el = sig13_el + f1_el_x*(z2 - z1);
- sig23_el = sig23_el + f1_el_y*(z2 - z1);
-
- }
-
- else
-
{
- Vector3r l =
- std::min(currentContactGeometry->radius2,
- currentContactGeometry->radius1)
- *currentContactGeometry->normal;
-
- /// Calcul des contraintes elastiques spheres/parois
-
- sig11_el = sig11_el + f1_el_x*l[0];
- sig22_el = sig22_el + f1_el_y*l[1];
- sig33_el = sig33_el + f1_el_z*l[2];
- sig12_el = sig12_el + f1_el_x*l[1];
- sig13_el = sig13_el + f1_el_x*l[2];
- sig23_el = sig23_el + f1_el_y*l[2];
-
- }
+ j=j+1;
- /// fabric tensor
+ Vector3r fel = currentContactPhysics->normalForce + currentContactPhysics->shearForce;
- Vector3r normal = currentContactGeometry->normal;
-
- for (int i=0; i<3; ++i)
- {
- for (int n=0; n<3; ++n)
- {
- //fabricTensor[i][n]
- FT[i][n]
- += normal[i]*normal[n];
+ f1_el_x=fel[0];
+ f1_el_y=fel[1];
+ f1_el_z=fel[2];
+
+ int geometryIndex1 =
+ ( *bodies ) [id1]->geometricalModel->getClassIndex();
+ int geometryIndex2 =
+ ( *bodies ) [id2]->geometricalModel->getClassIndex();
+
+ if ( geometryIndex1 == geometryIndex2 )
+
+ {
+ BodyMacroParameters* de1 = static_cast<BodyMacroParameters*> ( ( *bodies ) [id1]->physicalParameters.get() );
+ x1 = de1->se3.position[0];
+ y1 = de1->se3.position[1];
+ z1 = de1->se3.position[2];
+
+
+ BodyMacroParameters* de2 = static_cast<BodyMacroParameters*> ( ( *bodies ) [id2]->physicalParameters.get() );
+ x2 = de2->se3.position[0];
+ y2 = de2->se3.position[1];
+ z2 = de2->se3.position[2];
+
+ ///Calcul des contraintes elastiques spheres/spheres
+
+ sig11_el = sig11_el + f1_el_x* ( x2 - x1 );
+ sig22_el = sig22_el + f1_el_y* ( y2 - y1 );
+ sig33_el = sig33_el + f1_el_z* ( z2 - z1 );
+ sig12_el = sig12_el + f1_el_x* ( y2 - y1 );
+ sig13_el = sig13_el + f1_el_x* ( z2 - z1 );
+ sig23_el = sig23_el + f1_el_y* ( z2 - z1 );
+
}
+
+ else
+
+ {
+ Vector3r l = std::min ( currentContactGeometry->radius2, currentContactGeometry->radius1 ) *currentContactGeometry->normal;
+
+ /// Calcul des contraintes elastiques spheres/parois
+
+ sig11_el = sig11_el + f1_el_x*l[0];
+ sig22_el = sig22_el + f1_el_y*l[1];
+ sig33_el = sig33_el + f1_el_z*l[2];
+ sig12_el = sig12_el + f1_el_x*l[1];
+ sig13_el = sig13_el + f1_el_x*l[2];
+ sig23_el = sig23_el + f1_el_y*l[2];
+
+ }
+
+ /// fabric tensor
+
+ Vector3r normal = currentContactGeometry->normal;
+
+ for ( int i=0; i<3; ++i )
+ {
+ for ( int n=0; n<3; ++n )
+ {
+ //fabricTensor[i][n]
+ FT[i][n]
+ += normal[i]*normal[n];
+ }
+ }
+
}
-
- }
- }
- }
+ }
+ }
/// FabricTensor
-
+
//Real traceFT = (FT[0][0]+FT[1][1]+FT[2][2])/j;
-
+
/// calcul de l'energie cinetique:
Real nbElt = 0, SR = 0, Vs=0, Rbody=0, Rmin=1, Rmax=0;
-
+
BodyContainer::iterator bi = bodies->begin();
BodyContainer::iterator biEnd = bodies->end();
-
- for ( ; bi!=biEnd; ++bi)
-
- {
+
+ for ( ; bi!=biEnd; ++bi )
+
+ {
shared_ptr<Body> b = *bi;
-
+
int geometryIndex = b->geometricalModel->getClassIndex();
-
- if (geometryIndex == SpheresClassIndex)
+
+ if ( geometryIndex == SpheresClassIndex )
{
nbElt +=1;
- ParticleParameters* pp =
- static_cast<ParticleParameters*>(b->physicalParameters.get());
+ ParticleParameters* pp = static_cast<ParticleParameters*> ( b->physicalParameters.get() );
Vector3r v = pp->velocity;
kinematicE +=
- 0.5*(pp->mass)*(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
-
- Sphere* sphere =
- static_cast<Sphere*>(b->geometricalModel.get());
+ 0.5* ( pp->mass ) * ( v[0]*v[0]+v[1]*v[1]+v[2]*v[2] );
+
+ Sphere* sphere = static_cast<Sphere*> ( b->geometricalModel.get() );
Rbody = sphere->radius;
- if (Rbody<Rmin) Rmin = Rbody;
- if (Rbody>Rmax) Rmax = Rbody;
+ if ( Rbody<Rmin ) Rmin = Rbody;
+ if ( Rbody>Rmax ) Rmax = Rbody;
SR+=Rbody;
-
- Vs += 1.3333333*Mathr::PI*(Rbody*Rbody*Rbody);
-
+
+ Vs += 1.3333333*Mathr::PI* ( Rbody*Rbody*Rbody );
+
}
}
/// coordination number
-
+
Real coordN = 0;
- coordN = 2*(j/nbElt); // ????????????????????????????????????????????
-
+ coordN = 2* ( j/nbElt ); // ????????????????????????????????????????????
+
/// Calcul des contraintes "globales"
-
+
Real SIG_11_el=0, SIG_22_el=0, SIG_33_el=0, SIG_12_el=0, SIG_13_el=0, SIG_23_el=0;
-
+
// volume de l'echantillon
-
-// PhysicalParameters* p_bottom = static_cast<PhysicalParameters*>((*bodies)[wall_bottom_id]->physicalParameters.get());
-// PhysicalParameters* p_top = static_cast<PhysicalParameters*>((*bodies)[wall_top_id]->physicalParameters.get());
-// PhysicalParameters* p_left = static_cast<PhysicalParameters*>((*bodies)[wall_left_id]->physicalParameters.get());
-// PhysicalParameters* p_right = static_cast<PhysicalParameters*>((*bodies)[wall_right_id]->physicalParameters.get());
-// PhysicalParameters* p_front = static_cast<PhysicalParameters*>((*bodies)[wall_front_id]->physicalParameters.get());
-// PhysicalParameters* p_back = static_cast<PhysicalParameters*>((*bodies)[wall_back_id]->physicalParameters.get());
-
-
-// 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;
-
+
+// PhysicalParameters* p_bottom = static_cast<PhysicalParameters*>((*bodies)[wall_bottom_id]->physicalParameters.get());
+// PhysicalParameters* p_top = static_cast<PhysicalParameters*>((*bodies)[wall_top_id]->physicalParameters.get());
+// PhysicalParameters* p_left = static_cast<PhysicalParameters*>((*bodies)[wall_left_id]->physicalParameters.get());
+// PhysicalParameters* p_right = static_cast<PhysicalParameters*>((*bodies)[wall_right_id]->physicalParameters.get());
+// PhysicalParameters* p_front = static_cast<PhysicalParameters*>((*bodies)[wall_front_id]->physicalParameters.get());
+// PhysicalParameters* p_back = static_cast<PhysicalParameters*>((*bodies)[wall_back_id]->physicalParameters.get());
+
+
+// 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;
+
Real Rmoy = SR/nbElt;
- if (Omega::instance().getCurrentIteration() == 1)
- {cerr << "Rmoy = "<< Rmoy << " Rmin = " << Rmin << " Rmax = " <<
- Rmax << endl;}
-
+ if ( Omega::instance().getCurrentIteration() == 1 )
+ {
+ cerr << "Rmoy = "<< Rmoy << " Rmin = " << Rmin << " Rmax = " <<
+ Rmax << endl;
+ }
+
Real V = ( triaxCompEng->height ) * ( triaxCompEng->width ) * ( triaxCompEng->depth );
-
+
SIG_11_el = sig11_el/V;
SIG_22_el = sig22_el/V;
SIG_33_el = sig33_el/V;
SIG_12_el = sig12_el/V;
SIG_13_el = sig13_el/V;
SIG_23_el = sig23_el/V;
-
+
/// calcul des deformations
-
-// Real EPS_11=0, EPS_22=0, EPS_33=0;
-//
-// Real width_0 = upperCorner[0]-lowerCorner[0], height_0 =
-// upperCorner[1]-lowerCorner[1],
-// depth_0 = upperCorner[2]-lowerCorner[2];
-//
-// // EPS_11 = (width_0 - width)/width_0;
-// EPS_11 = std::log(width_0) - std::log(width);
-// // EPS_22 = (height_0 - height)/height_0;
-// EPS_22 = std::log(height_0) - std::log(height);
-// // EPS_33 = (depth_0 - depth)/depth_0;
-// EPS_33 = std::log(depth_0) - std::log(depth);
-
- /// porosity
-
+
+// Real EPS_11=0, EPS_22=0, EPS_33=0;
+//
+// Real width_0 = upperCorner[0]-lowerCorner[0], height_0 =
+// upperCorner[1]-lowerCorner[1],
+// depth_0 = upperCorner[2]-lowerCorner[2];
+//
+// // EPS_11 = (width_0 - width)/width_0;
+// EPS_11 = std::log(width_0) - std::log(width);
+// // EPS_22 = (height_0 - height)/height_0;
+// EPS_22 = std::log(height_0) - std::log(height);
+// // EPS_33 = (depth_0 - depth)/depth_0;
+// EPS_33 = std::log(depth_0) - std::log(depth);
+
+ /// porosity
+
Real Vv = V - Vs;
-
+
Real n = Vv/V;
-// Real e = Vv/Vs;
+// Real e = Vv/Vs;
-// mise a zero des deformations qd comp triaxiale commence
-// if (triaxCompEng->compressionActivated) { }
-
+// mise a zero des deformations qd comp triaxiale commence
+// if (triaxCompEng->compressionActivated) { }
+
/// UnbalancedForce
-
- Real equilibriumForce = triaxCompEng->ComputeUnbalancedForce(ncb);
-// Real equilibriumForce = sampleCapPressEng->ComputeUnbalancedForce(body);
- if (Omega::instance().getCurrentIteration() % 100 == 0)
- {cerr << "current Iteration " << Omega::instance().getCurrentIteration()
- << endl;}
-
- ofile /*<< lexical_cast<string>(Omega::instance().getSimulationTime()) << " "*/
- << lexical_cast<string>(Omega::instance().getCurrentIteration()) << " "
- << lexical_cast<string>(kinematicE)<< " "
- << lexical_cast<string>(equilibriumForce)<<" "
- << lexical_cast<string>(n)<<" "
-// << lexical_cast<string>(e)<<" "
- << lexical_cast<string>(coordN)<<" "
- << lexical_cast<string>(SIG_11_el) << " "
- << lexical_cast<string>(SIG_22_el) << " "
- << lexical_cast<string>(SIG_33_el) << " "
- << lexical_cast<string>(SIG_12_el) << " "
- << lexical_cast<string>(SIG_13_el)<< " "
- << lexical_cast<string>(SIG_23_el)<< " "
- << lexical_cast<string>(triaxCompEng->strain[0]) << " "
- << lexical_cast<string>(triaxCompEng->strain[1]) << " "
- << lexical_cast<string>(triaxCompEng->strain[2]) << " "
- << lexical_cast<string>(FT[0][0]/j)<<" "
- << lexical_cast<string>(FT[0][1]/j)<<" "
- << lexical_cast<string>(FT[0][2]/j)<<" "
- << lexical_cast<string>(FT[1][0]/j)<<" "
- << lexical_cast<string>(FT[1][1]/j)<<" "
- << lexical_cast<string>(FT[1][2]/j)<<" "
- << lexical_cast<string>(FT[2][0]/j)<<" "
- << lexical_cast<string>(FT[2][1]/j)<<" "
- << lexical_cast<string>(FT[2][2]/j)<<" "
+ Real equilibriumForce = triaxCompEng->ComputeUnbalancedForce ( ncb );
+// Real equilibriumForce = sampleCapPressEng->ComputeUnbalancedForce(body);
+
+ if ( Omega::instance().getCurrentIteration() % 100 == 0 )
+ {
+ cerr << "current Iteration " << Omega::instance().getCurrentIteration()
<< endl;
-
+ }
+
+ ofile /*<< lexical_cast<string>(Omega::instance().getSimulationTime()) << " "*/
+ << lexical_cast<string> ( Omega::instance().getCurrentIteration() ) << " "
+ << lexical_cast<string> ( kinematicE ) << " "
+ << lexical_cast<string> ( equilibriumForce ) <<" "
+ << lexical_cast<string> ( n ) <<" "
+// << lexical_cast<string>(e)<<" "
+ << lexical_cast<string> ( coordN ) <<" "
+ << lexical_cast<string> ( SIG_11_el ) << " "
+ << lexical_cast<string> ( SIG_22_el ) << " "
+ << lexical_cast<string> ( SIG_33_el ) << " "
+ << lexical_cast<string> ( SIG_12_el ) << " "
+ << lexical_cast<string> ( SIG_13_el ) << " "
+ << lexical_cast<string> ( SIG_23_el ) << " "
+ << lexical_cast<string> ( triaxCompEng->strain[0] ) << " "
+ << lexical_cast<string> ( triaxCompEng->strain[1] ) << " "
+ << lexical_cast<string> ( triaxCompEng->strain[2] ) << " "
+ << lexical_cast<string> ( FT[0][0]/j ) <<" "
+ << lexical_cast<string> ( FT[0][1]/j ) <<" "
+ << lexical_cast<string> ( FT[0][2]/j ) <<" "
+ << lexical_cast<string> ( FT[1][0]/j ) <<" "
+ << lexical_cast<string> ( FT[1][1]/j ) <<" "
+ << lexical_cast<string> ( FT[1][2]/j ) <<" "
+ << lexical_cast<string> ( FT[2][0]/j ) <<" "
+ << lexical_cast<string> ( FT[2][1]/j ) <<" "
+ << lexical_cast<string> ( FT[2][2]/j ) <<" "
+ << endl;
+
}
YADE_PLUGIN();
Added: trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.cpp 2009-06-24 16:28:58 UTC (rev 1812)
+++ trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.cpp 2009-06-24 17:04:39 UTC (rev 1813)
@@ -0,0 +1,152 @@
+/*************************************************************************
+* Copyright (C) 2006 by Luc Scholtes *
+* luc.scholtes@xxxxxxxxxxx *
+* *
+* 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 "SampleCapillaryPressureEngine.hpp"
+#include <yade/pkg-dem/CapillaryCohesiveLaw.hpp>
+#include<yade/core/MetaBody.hpp>
+#include<yade/core/Omega.hpp>
+//#include<yade/pkg-common/Force.hpp>
+#include<yade/pkg-dem/ElasticContactInteraction.hpp>
+#include<yade/lib-base/yadeWm3Extra.hpp>
+#include <boost/lexical_cast.hpp>
+
+using namespace boost;
+using namespace std;
+
+SampleCapillaryPressureEngine::SampleCapillaryPressureEngine()
+{
+ //cerr << "constructeur de SamplePressureEngine" << endl;
+
+ capillaryCohesiveLaw = new CapillaryCohesiveLaw;
+ capillaryCohesiveLaw->CapillaryPressure= 0;
+ capillaryCohesiveLaw->sdecGroupMask = 2;
+
+ StabilityCriterion=0.01;
+ SigmaPrecision = 0.001;
+ Phase1=false;
+ Phase1End = "Phase1End";
+ Iteration = 0;
+ UnbalancedForce = 0.01;
+
+ pressureVariationActivated=false;
+ Pressure = 0;
+ PressureVariation = 1000;
+ fusionDetection = true;
+ binaryFusion = true;
+
+ pressureIntervalRec = 10000;
+
+ //cerr << "fin du constructeur de SamplePressureEngine" << endl;
+
+}
+
+SampleCapillaryPressureEngine::~SampleCapillaryPressureEngine()
+{
+}
+
+void SampleCapillaryPressureEngine::registerAttributes()
+{
+ TriaxialStressController::registerAttributes();
+ REGISTER_ATTRIBUTE(PressureVariation);
+ REGISTER_ATTRIBUTE(Pressure);
+ REGISTER_ATTRIBUTE(UnbalancedForce);
+ REGISTER_ATTRIBUTE(StabilityCriterion);
+ REGISTER_ATTRIBUTE(SigmaPrecision);
+ REGISTER_ATTRIBUTE(pressureVariationActivated);
+ REGISTER_ATTRIBUTE(fusionDetection);
+ REGISTER_ATTRIBUTE(binaryFusion);
+ REGISTER_ATTRIBUTE(pressureIntervalRec);
+
+}
+
+void SampleCapillaryPressureEngine::updateParameters(MetaBody * ncb)
+{
+ //cerr << "updateParameters" << endl;
+// MetaBody * ncb = static_cast<MetaBody*>(body);
+
+ UnbalancedForce = ComputeUnbalancedForce(ncb);
+
+ if (Omega::instance().getCurrentIteration() % 100 == 0) cerr << "UnbalancedForce=" << UnbalancedForce << endl;
+
+ if (!Phase1 && UnbalancedForce<=StabilityCriterion && !pressureVariationActivated)
+ {
+ internalCompaction = false;
+ Phase1 = true;
+ }
+
+ if ( Phase1 && UnbalancedForce<=StabilityCriterion && !pressureVariationActivated)
+
+ {
+ Real S = meanStress; // revoir ici comment acceder à computeStress(ncb);
+ //abs ( (meanStress-sigma_iso ) /sigma_iso ) <0.005
+
+ cerr << "Smoy = " << meanStress << endl;
+ if ((S > (sigma_iso - (sigma_iso*SigmaPrecision))) && (S < (sigma_iso + (sigma_iso*SigmaPrecision))))
+ {
+ //Iteration = Omega::instance().getCurrentIteration();
+
+ // saving snapshot.xml
+ string fileName = "../data/" + Phase1End + "_" +
+ lexical_cast<string>(Omega::instance().getCurrentIteration()) + ".xml";
+ cerr << "saving snapshot: " << fileName << " ...";
+ Omega::instance().saveSimulation(fileName);
+
+ //recordStructure(ncb, Omega::instance().getCurrentIteration());
+
+ pressureVariationActivated = true;
+ }
+
+ }
+
+// if (pressureVariationActivated)
+//
+// {
+// autoCompressionActivation = true;
+// }
+
+ //cerr << "fin updateParameters" << endl;
+}
+
+void SampleCapillaryPressureEngine::applyCondition(MetaBody * ncb)
+{
+ //cerr << "applyConditionSampleCapillaryPressure" << endl;
+// MetaBody* ncb = static_cast<MetaBody*>(body);
+
+ updateParameters(ncb);
+
+ TriaxialStressController::applyCondition(ncb);
+
+ //cerr << "1" << endl;
+
+ if (pressureVariationActivated)
+
+ {
+ //if ((Omega::instance().getCurrentIteration()) % pressureIntervalRec == 0) recordStructure(ncb, Omega::instance().getCurrentIteration());
+
+ if (Omega::instance().getCurrentIteration() % 100 == 0) cerr << "pressure variation!!" << endl;
+
+ if ((Pressure>=0) && (Pressure<=1000000000)) Pressure += PressureVariation;
+ capillaryCohesiveLaw->CapillaryPressure = Pressure;
+
+ capillaryCohesiveLaw->fusionDetection = fusionDetection;
+ capillaryCohesiveLaw->binaryFusion = binaryFusion;
+
+ }
+
+ else { capillaryCohesiveLaw->CapillaryPressure = Pressure;
+ capillaryCohesiveLaw->fusionDetection = fusionDetection;
+ capillaryCohesiveLaw->binaryFusion = binaryFusion;}
+
+ if (Omega::instance().getCurrentIteration() % 100 == 0) cerr << "capillary pressure = " << Pressure << endl;
+
+ //cerr << "3" << endl;
+ capillaryCohesiveLaw->action(ncb);
+ //cerr << "4" << endl;
+
+ UnbalancedForce = ComputeUnbalancedForce(ncb);
+}
Property changes on: trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.cpp
___________________________________________________________________
Name: svn:executable
+ *
Added: trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.hpp 2009-06-24 16:28:58 UTC (rev 1812)
+++ trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.hpp 2009-06-24 17:04:39 UTC (rev 1813)
@@ -0,0 +1,65 @@
+/*************************************************************************
+* Copyright (C) 2006 by Luc Scholtes *
+* luc.scholtes@xxxxxxxxxxx *
+* *
+* 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 SAMPLE_CAPILLARY_PRESSURE_ENGINE_HPP
+#define SAMPLE_CAPILLARY_PRESSURE_ENGINE_HPP
+
+#include<yade/core/DeusExMachina.hpp>
+#include <Wm3Vector3.h>
+#include<yade/lib-base/yadeWm3.hpp>
+#include "TriaxialStressController.hpp"
+#include <string>
+
+/*! \brief Isotropic compression + suction variation test */
+
+class CapillaryCohesiveLaw;
+class PhysicalAction;
+
+
+class SampleCapillaryPressureEngine : public TriaxialStressController
+{
+ private :
+ //shared_ptr<PhysicalAction> actionForce;
+
+ public :
+ SampleCapillaryPressureEngine();
+ virtual ~SampleCapillaryPressureEngine();
+
+ unsigned int interval, VariationInterval;
+
+ //! Max ratio of resultant forces on mean contact force
+ Real UnbalancedForce;
+ //! Value of UnbalancedForce for which the system is considered
+ Real StabilityCriterion, SigmaPrecision;
+ //! is isotropicInternalCompactionFinished?
+ bool Phase1;
+ int Iteration, pressureIntervalRec;
+ std::string Phase1End;
+ //! pressure affectation
+ Real Pressure;
+ Real PressureVariation;
+ //! Is pressure variation currently activated?
+ bool pressureVariationActivated, fusionDetection, binaryFusion;
+
+ //shared_ptr<CapillaryCohesiveLaw> capillaryCohesiveLaw;
+ CapillaryCohesiveLaw* capillaryCohesiveLaw;
+
+ virtual void applyCondition(MetaBody * ncb);
+ void updateParameters(MetaBody * ncb);
+
+
+ protected :
+ virtual void registerAttributes();
+ REGISTER_CLASS_NAME(SampleCapillaryPressureEngine);
+ REGISTER_BASE_CLASS_NAME(TriaxialStressController);
+};
+
+REGISTER_SERIALIZABLE(SampleCapillaryPressureEngine);
+
+#endif // SAMPLECAPILLARYPRESSUREENGINE_HPP
+
Property changes on: trunk/pkg/dem/Engine/DeusExMachina/SampleCapillaryPressureEngine.hpp
___________________________________________________________________
Name: svn:executable
+ *
Modified: trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp 2009-06-24 16:28:58 UTC (rev 1812)
+++ trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp 2009-06-24 17:04:39 UTC (rev 1813)
@@ -9,6 +9,7 @@
//Modifs : Parameters renamed as MeniscusParameters
//id1/id2 classés pour que id1 soit toujours le plus petit grain, FIXME : angle de mouillage?
//FIXME : dans triaxialStressController, changer le test de nullité de la force dans updateStiffness
+//FIXME : needs "requestErase" somewhere
#include "CapillaryCohesiveLaw.hpp"
Modified: trunk/pkg/dem/PreProcessor/TriaxialTestWater.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TriaxialTestWater.cpp 2009-06-24 16:28:58 UTC (rev 1812)
+++ trunk/pkg/dem/PreProcessor/TriaxialTestWater.cpp 2009-06-24 17:04:39 UTC (rev 1813)
@@ -18,7 +18,7 @@
#include<yade/pkg-dem/ElasticContactLaw.hpp>
#include <yade/pkg-dem/CapillaryCohesiveLaw.hpp>
// #include<yade/pkg-dem/SimpleElasticRelationships.hpp>
-/////////////#include<yade/pkg-dem/SimpleElasticRelationshipsWater.hpp>
+#include<yade/pkg-dem/SimpleElasticRelationshipsWater.hpp>
#include<yade/pkg-dem/BodyMacroParameters.hpp>
#include<yade/pkg-dem/SDECLinkPhysics.hpp>
@@ -44,6 +44,7 @@
#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
#include<yade/pkg-common/MetaInteractingGeometry2AABB.hpp>
#include<yade/pkg-common/MetaInteractingGeometry.hpp>
+#include<yade/pkg-common/InteractingSphere2AABB.hpp>
#include<yade/pkg-common/GravityEngines.hpp>
#include<yade/pkg-common/PhysicalActionApplier.hpp>
@@ -57,6 +58,7 @@
#include<yade/core/Body.hpp>
#include<yade/pkg-common/InteractingBox.hpp>
#include<yade/pkg-common/InteractingSphere.hpp>
+#include<yade/pkg-dem/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp>
#include<yade/pkg-common/PhysicalActionContainerReseter.hpp>
@@ -503,9 +505,17 @@
void TriaxialTestWater::createActors(shared_ptr<MetaBody>& rootBody)
{
+ Real distanceFactor = 1.3;//Create potential interactions as soon as the distance is less than factor*(rad1+rad2)
+
shared_ptr<InteractionGeometryMetaEngine> interactionGeometryDispatcher(new InteractionGeometryMetaEngine);
- interactionGeometryDispatcher->add("InteractingSphere2InteractingSphere4SpheresContactGeometry");
+
+ shared_ptr<InteractingSphere2InteractingSphere4SpheresContactGeometry> iS2IS4SContactGeometry(new InteractingSphere2InteractingSphere4SpheresContactGeometry);
+ iS2IS4SContactGeometry->interactionDetectionFactor = distanceFactor;//Detect potential distant interaction (meniscii)
+
+ interactionGeometryDispatcher->add(iS2IS4SContactGeometry);
interactionGeometryDispatcher->add("InteractingBox2InteractingSphere4SpheresContactGeometry");
+
+
shared_ptr<InteractionPhysicsMetaEngine> interactionPhysicsDispatcher(new InteractionPhysicsMetaEngine);
// interactionPhysicsDispatcher->add("SimpleElasticRelationships");
@@ -513,12 +523,16 @@
/// OLD
//interactionPhysicsDispatcher->add("BodyMacroParameters","BodyMacroParameters","MacroMicroElasticRelationshipsWater");
/// NEW
-///////////// shared_ptr<InteractionPhysicsEngineUnit> ss(new SimpleElasticRelationshipsWater);
-///////////// interactionPhysicsDispatcher->add(ss);
+ shared_ptr<InteractionPhysicsEngineUnit> ss(new SimpleElasticRelationshipsWater);
+ interactionPhysicsDispatcher->add(ss);
shared_ptr<BoundingVolumeMetaEngine> boundingVolumeDispatcher = shared_ptr<BoundingVolumeMetaEngine>(new BoundingVolumeMetaEngine);
- boundingVolumeDispatcher->add("InteractingSphere2AABB");
+
+ shared_ptr<InteractingSphere2AABB> interactingSphere2AABB(new InteractingSphere2AABB);
+ interactingSphere2AABB->aabbEnlargeFactor = distanceFactor;//Detect potential distant interaction (meniscii)
+
+ boundingVolumeDispatcher->add(interactingSphere2AABB);
boundingVolumeDispatcher->add("InteractingBox2AABB");
boundingVolumeDispatcher->add("MetaInteractingGeometry2AABB");
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-06-24 16:28:58 UTC (rev 1812)
+++ trunk/pkg/dem/SConscript 2009-06-24 17:04:39 UTC (rev 1813)
@@ -3,23 +3,23 @@
import os.path,os
# Dir('#') expands to the root of the source tree -- that is done by scons
-if os.path.exists(Dir('#').abspath+'/extra/triangulation/libTesselationWrapper.a'):
- print "Will build VolumetricContactLaw since libTesselationWrapper.a was found."
- env.Install('$PREFIX/lib/yade$SUFFIX/pkg-dem',[
- env.SharedLibrary('VolumicContactLaw',
- ['Engine/StandAloneEngine/VolumicContactLaw.cpp'],
- LIBPATH=env['LIBPATH']+['../extra/triangulation'],
- LIBS=env['LIBS']+['SDECLinkPhysics',
- 'ElasticContactInteraction',
- 'SDECLinkGeometry',
- 'SpheresContactGeometry',
- 'BodyMacroParameters',
- 'Sphere',
- 'RigidBodyParameters',
- 'InteractingSphere',
- 'TesselationWrapper',
- 'CGAL'])
- ])
+#if os.path.exists(Dir('#').abspath+'/extra/triangulation/libTesselationWrapper.a'):
+ #print "Will build VolumetricContactLaw since libTesselationWrapper.a was found."
+ #env.Install('$PREFIX/lib/yade$SUFFIX/pkg-dem',[
+ #env.SharedLibrary('VolumicContactLaw',
+ #['Engine/StandAloneEngine/VolumicContactLaw.cpp'],
+ #LIBPATH=env['LIBPATH']+['../extra/triangulation'],
+ #LIBS=env['LIBS']+['SDECLinkPhysics',
+ #'ElasticContactInteraction',
+ #'SDECLinkGeometry',
+ #'SpheresContactGeometry',
+ #'BodyMacroParameters',
+ #'Sphere',
+ #'RigidBodyParameters',
+ #'InteractingSphere',
+ #'TesselationWrapper',
+ #'CGAL'])
+ #])
cpmModel='../brefcom-mm.hh'
if os.path.exists('../../'+cpmModel):
@@ -313,7 +313,6 @@
'ElasticContactInteraction',
'SpheresContactGeometry',
- 'MacroMicroElasticRelationships',
'Sphere',
]),
@@ -793,7 +792,7 @@
TriaxialStateRecorder
CapillaryStressRecorder
ContactStressRecorder
- MacroMicroElasticRelationships
+ SimpleElasticRelationshipsWater
ElasticCriterionTimeStepper
InteractingSphere
InteractingBox
@@ -806,6 +805,7 @@
PhysicalActionContainerReseter
InteractionGeometryMetaEngine
InteractionPhysicsMetaEngine
+ InteractingSphere2InteractingSphere4SpheresContactGeometry
PhysicalActionApplier
PhysicalParametersMetaEngine
BoundingVolumeMetaEngine
@@ -815,8 +815,8 @@
AABB
PersistentSAPCollider
MetaInteractingGeometry2AABB
+ InteractingSphere2AABB
TriaxialStressController
- MacroMicroElasticRelationshipsWater
TriaxialCompressionEngine
GlobalStiffnessTimeStepper''')),
@@ -859,21 +859,21 @@
'InteractingSphere',
'InteractingBox' ]),
-env.SharedLibrary('MacroMicroElasticRelationshipsWater',
- ['Engine/EngineUnit/MacroMicroElasticRelationshipsWater.cpp'],
- LIBS=env['LIBS']+['SDECLinkPhysics',
- 'SDECLinkGeometry',
- 'ElasticContactInteraction',
- 'SpheresContactGeometry',
- 'BodyMacroParameters',
- 'RigidBodyParameters',
- 'ParticleParameters',
- 'InteractionPhysicsMetaEngine',
- 'CapillaryParameters']),
+#env.SharedLibrary('MacroMicroElasticRelationshipsWater',
+ #['Engine/EngineUnit/MacroMicroElasticRelationshipsWater.cpp'],
+ #LIBS=env['LIBS']+['SDECLinkPhysics',
+ #'SDECLinkGeometry',
+ #'ElasticContactInteraction',
+ #'SpheresContactGeometry',
+ #'BodyMacroParameters',
+ #'RigidBodyParameters',
+ #'ParticleParameters',
+ #'InteractionPhysicsMetaEngine',
+ #'CapillaryParameters']),
-env.SharedLibrary('InteractingSphere2AABBwater',
- ['Engine/EngineUnit/InteractingSphere2AABBwater.cpp'],
- LIBS=env['LIBS']+['InteractingSphere', 'AABB', 'BoundingVolumeMetaEngine']),
+#env.SharedLibrary('InteractingSphere2AABBwater',
+ #['Engine/EngineUnit/InteractingSphere2AABBwater.cpp'],
+ #LIBS=env['LIBS']+['InteractingSphere', 'AABB', 'BoundingVolumeMetaEngine']),
env.SharedLibrary('HydraulicForceEngine',
['Engine/DeusExMachina/HydraulicForceEngine.cpp'],
@@ -883,6 +883,15 @@
'$PREFIX/include',
'DataClass/PhysicalParameters']),
+env.SharedLibrary('SampleCapillaryPressureEngine',
+ ['Engine/DeusExMachina/SampleCapillaryPressureEngine.cpp'],
+ LIBS=env['LIBS']+['yade-factory',
+ 'yade-base',
+ 'CapillaryCohesiveLaw',
+ 'ParticleParameters',
+ 'ElasticContactInteraction',
+ 'TriaxialStressController']),
+
env.SharedLibrary('StaticSpheresAttractionEngine',
['Engine/DeusExMachina/StaticSpheresAttractionEngine.cpp'],
LIBS=env['LIBS']+['yade-base', 'Sphere', 'StaticAttractionEngine', 'ElasticContactInteraction'],
@@ -1079,9 +1088,21 @@
#'/home/bruno/YADE/trunk_svn/extra/triangulation/TriaxialState.cpp'],
#"make",
#chdir='/home/bruno/YADE/trunk_svn/extra/triangulation')
+
+ env.SharedLibrary('SimpleElasticRelationshipsWater',
+ ['Engine/EngineUnit/SimpleElasticRelationshipsWater.cpp'],
+ LIBS=env['LIBS']+['SDECLinkPhysics',
+ 'ElasticContactInteraction',
+ 'CapillaryParameters',
+ 'SpheresContactGeometry',
+ 'BodyMacroParameters',
+ 'RigidBodyParameters',
+ 'ParticleParameters',
+ 'InteractionPhysicsMetaEngine']),
])
+
if 'YADE_OPENGL' in env['CPPDEFINES']:
env.Install('$PREFIX/lib/yade$SUFFIX/pkg-dem',[
