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.schol...@hmg.inpg.fr * +* * +* 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.schol...@hmg.inpg.fr * +* * +* 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',[ _______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp