yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #08134
[Branch ~yade-dev/yade/trunk] Rev 2974: - capillary model: nullify Fcap and volume when P changes from something to zero
------------------------------------------------------------
revno: 2974
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Wed 2011-11-30 18:39:33 +0100
message:
- capillary model: nullify Fcap and volume when P changes from something to zero
- remove the if(!b) test in a few important places (proof of concept + don't waste time doing the same test twice)
modified:
pkg/common/GravityEngines.cpp
pkg/common/InteractionLoop.cpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp
pkg/dem/NewtonIntegrator.cpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/common/GravityEngines.cpp'
--- pkg/common/GravityEngines.cpp 2011-03-01 10:31:54 +0000
+++ pkg/common/GravityEngines.cpp 2011-11-30 17:39:33 +0000
@@ -19,7 +19,7 @@
const Real dt(scene->dt);
YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
// skip clumps, only apply forces on their constituents
- if(!b || b->isClump()) continue;
+ if(b->isClump()) continue;
if(mask!=0 && (b->groupMask & mask)==0) continue;
scene->forces.addForce(b->getId(),gravity*b->state->mass);
// work done by gravity is "negative", since the energy appears in the system from outside
@@ -30,7 +30,7 @@
void CentralGravityEngine::action(){
const Vector3r& centralPos=Body::byId(centralBody)->state->pos;
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if(!b || b->isClump() || b->getId()==centralBody) continue; // skip clumps and central body
+ if(b->isClump() || b->getId()==centralBody) continue; // skip clumps and central body
if(mask!=0 && (b->groupMask & mask)==0) continue;
Real F=accel*b->state->mass;
Vector3r toCenter=centralPos-b->state->pos; toCenter.normalize();
=== modified file 'pkg/common/InteractionLoop.cpp'
--- pkg/common/InteractionLoop.cpp 2011-02-27 14:26:15 +0000
+++ pkg/common/InteractionLoop.cpp 2011-11-30 17:39:33 +0000
@@ -72,7 +72,7 @@
#ifdef YADE_SUBDOMAINS
- YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){ if(unlikely(!b)) continue; FOREACH(const Body::MapId2IntrT::value_type& mapItem, b->intrs){
+ YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){ FOREACH(const Body::MapId2IntrT::value_type& mapItem, b->intrs){
const shared_ptr<Interaction>& I(mapItem.second);
#else
#ifdef YADE_OPENMP
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2011-11-30 13:30:02 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2011-11-30 17:39:33 +0000
@@ -31,7 +31,7 @@
void Law2_ScGeom_CapillaryPhys_Capillarity::postLoad(Law2_ScGeom_CapillaryPhys_Capillarity&){
- capillary = shared_ptr<capillarylaw>(new capillarylaw); // ????????
+ capillary = shared_ptr<capillarylaw>(new capillarylaw);
capillary->fill("M(r=1)");
capillary->fill("M(r=1.1)");
capillary->fill("M(r=1.25)");
@@ -67,93 +67,84 @@
void Law2_ScGeom_CapillaryPhys_Capillarity::action()
{
if (!scene) cerr << "scene not defined!";
- shared_ptr<BodyContainer>& bodies = scene->bodies;
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
+ if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene);
- if (fusionDetection) {
- if (!bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene);
- //bodiesMenisciiList.display();
- }
-
- InteractionContainer::iterator ii = scene->interactions->begin();
- InteractionContainer::iterator iiEnd = scene->interactions->end();
+ InteractionContainer::iterator ii = scene->interactions->begin();
+ InteractionContainer::iterator iiEnd = scene->interactions->end();
bool hertzInitialized = false;
- for( ; ii!=iiEnd ; ++ii ) {
- if ((*ii)->isReal()) {
- const shared_ptr<Interaction>& interaction = *ii;
+ for (; ii!=iiEnd ; ++ii) {
+ if ((*ii)->isReal()) {
+ const shared_ptr<Interaction>& interaction = *ii;
if (!hertzInitialized) {//NOTE: We are assuming that only one type is used in one simulation here
if (CapillaryPhys::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=false;
else if (MindlinCapillaryPhys::getClassIndexStatic()==interaction->phys->getClassIndex()) hertzOn=true;
- else LOG_ERROR("The capillary law is not implemented for interactions using"<<interaction->phys->getClassName());}
+ else LOG_ERROR("The capillary law is not implemented for interactions using"<<interaction->phys->getClassName());
+ }
hertzInitialized = true;
- unsigned int id1 = interaction->getId1();
- unsigned int id2 = interaction->getId2();
-
- Body* b1 = (*bodies)[id1].get();
- Body* b2 = (*bodies)[id2].get();
-
- if (b1 and b2) {// FIXME: possible memory leak here?
-
- /// interaction geometry search (this test is to compute capillarity only between spheres (probably a better way to do that)
- int geometryIndex1 = (*bodies)[id1]->shape->getClassIndex(); // !!!
- int geometryIndex2 = (*bodies)[id2]->shape->getClassIndex();
- if (!(geometryIndex1 == geometryIndex2)) continue;
-
- /// definition of interacting objects (not necessarily in contact)
- ScGeom* currentContactGeometry = static_cast<ScGeom*>(interaction->geom.get());
- //CapillaryPhys* currentContactPhysics = static_cast<CapillaryPhys*>(interaction->phys.get());//obsolete
-
- /// contact physics depends on the contact law, that is used (either linear model or hertz model)
- CapillaryPhys* cundallContactPhysics=NULL;
- MindlinCapillaryPhys* mindlinContactPhysics=NULL;
- if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys*>(interaction->phys.get());//use CapillaryPhys for linear model
- else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>(interaction->phys.get());//use MindlinCapillaryPhys for hertz model
-
- /// Capillary components definition:
- Real liquidTension = 0.073; // superficial water tension N/m (0.073 is water tension at 20 Celsius degrees)
- //Real teta = 0; // perfect wetting (as in the case of pure water and glass beads)
-
- /// Interacting Grains:
- // If you want to define a ratio between YADE sphere size and real sphere size (Rk: OK if no gravity is accounted for)
- Real alpha=1;
- Real R1 = 0;
- R1=alpha*std::min(currentContactGeometry->radius2,currentContactGeometry->radius1 ) ;
- Real R2 = 0;
- R2=alpha*std::max(currentContactGeometry->radius2,currentContactGeometry->radius1 ) ;
- //cerr << "R1 = " << R1 << " R2 = "<< R2 << endl;
-
- /// intergranular distance
- Real D = alpha*((b2->state->pos-b1->state->pos).norm()-(currentContactGeometry->radius1+ currentContactGeometry->radius2)); // scGeom->penetrationDepth could probably be used here?
-
- if ((currentContactGeometry->penetrationDepth>=0)|| D<=0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere
- D=0; // defines Fcap when spheres interpenetrate //FIXME : D<0 leads to wrong interpolation has D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
- if (!hertzOn){
- if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
- cundallContactPhysics->meniscus=true;
- } else {
- if (fusionDetection && !mindlinContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
- mindlinContactPhysics->meniscus=true;}
- }
- Real Dinterpol = D/R2;
-
- /// Suction (Capillary pressure):
- Real Pinterpol =CapillaryPressure*(R2/liquidTension);
- if (!hertzOn) cundallContactPhysics->CapillaryPressure = CapillaryPressure;
- else mindlinContactPhysics->CapillaryPressure = CapillaryPressure;
-
- /// Capillary solution finder:
-// if (!hertzOn) {
+ unsigned int id1 = interaction->getId1();
+ unsigned int id2 = interaction->getId2();
+ Body* b1 = (*bodies)[id1].get();
+ Body* b2 = (*bodies)[id2].get();
+
+// /*if (b1 and b2)*/ {// FIXME: possible memory leak here?
+ /// interaction geometry search (this test is to compute capillarity only between spheres (probably a better way to do that)
+ int geometryIndex1 = (*bodies)[id1]->shape->getClassIndex(); // !!!
+ int geometryIndex2 = (*bodies)[id2]->shape->getClassIndex();
+ if (!(geometryIndex1 == geometryIndex2)) continue;
+
+ /// definition of interacting objects (not necessarily in contact)
+ ScGeom* currentContactGeometry = static_cast<ScGeom*>(interaction->geom.get());
+
+ /// contact physics depends on the contact law, that is used (either linear model or hertz model)
+ CapillaryPhys* cundallContactPhysics=NULL;
+ MindlinCapillaryPhys* mindlinContactPhysics=NULL;
+ if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys*>(interaction->phys.get());//use CapillaryPhys for linear model
+ else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>(interaction->phys.get());//use MindlinCapillaryPhys for hertz model
+
+ /// Capillary components definition:
+ Real liquidTension = 0.073; // superficial water tension N/m (0.073 is water tension at 20 Celsius degrees)
+
+ /// Interacting Grains:
+ // If you want to define a ratio between YADE sphere size and real sphere size
+ Real alpha=1;
+ Real R1 = 0;
+ R1=alpha*std::min(currentContactGeometry->radius2,currentContactGeometry->radius1) ;
+ Real R2 = 0;
+ R2=alpha*std::max(currentContactGeometry->radius2,currentContactGeometry->radius1) ;
+
+ /// intergranular distance
+ Real D = alpha*((b2->state->pos-b1->state->pos).norm()-(currentContactGeometry->radius1+ currentContactGeometry->radius2)); // scGeom->penetrationDepth could probably be used here?
+
+ if ((currentContactGeometry->penetrationDepth>=0)|| D<=0 || createDistantMeniscii) { //||(scene->iter < 1) ) // a simplified way to define meniscii everywhere
+ D=0; // defines Fcap when spheres interpenetrate. D<0 leads to wrong interpolation has D<0 has no solution in the interpolation : this is not physically interpretable!! even if, interpenetration << grain radius.
+ if (!hertzOn) {
+ if (fusionDetection && !cundallContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
+ cundallContactPhysics->meniscus=true;
+ } else {
+ if (fusionDetection && !mindlinContactPhysics->meniscus) bodiesMenisciiList.insert((*ii));
+ mindlinContactPhysics->meniscus=true;
+ }
+ }
+ Real Dinterpol = D/R2;
+
+ /// Suction (Capillary pressure):
+ Real Pinterpol =CapillaryPressure*(R2/liquidTension);
+ if (!hertzOn) cundallContactPhysics->CapillaryPressure = CapillaryPressure;
+ else mindlinContactPhysics->CapillaryPressure = CapillaryPressure;
+
+ /// Capillary solution finder:
if ((Pinterpol>=0) && (hertzOn? mindlinContactPhysics->meniscus : cundallContactPhysics->meniscus)) {
int* currentIndexes = hertzOn? mindlinContactPhysics->currentIndexes : cundallContactPhysics->currentIndexes;
+ //If P=0, we use null solution
MeniscusParameters
- solution(capillary->Interpolate(R1,R2,Dinterpol, Pinterpol, currentIndexes));
-
+ solution(Pinterpol? capillary->Interpolate(R1,R2,Dinterpol, Pinterpol, currentIndexes) : MeniscusParameters());
/// capillary adhesion force
Real Finterpol = solution.F;
Vector3r Fcap = Finterpol*(2*Mathr::PI*(R2/alpha)*liquidTension)*currentContactGeometry->normal;
if (!hertzOn) cundallContactPhysics->Fcap = Fcap;
- else mindlinContactPhysics->Fcap = Fcap;
-
+ else mindlinContactPhysics->Fcap = Fcap;
/// meniscus volume
Real Vinterpol = solution.V * pow(R2/alpha,3);
if (!hertzOn) {
@@ -165,27 +156,25 @@
if (Vinterpol != 0) mindlinContactPhysics->meniscus = true;
else mindlinContactPhysics->meniscus = false;
}
- if (!Vinterpol){
+ if (!Vinterpol) {
if (fusionDetection) bodiesMenisciiList.remove((*ii));
if (D>0) scene->interactions->requestErase(id1,id2);
}
-
/// wetting angles
if (!hertzOn) {
- cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
+ cundallContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
cundallContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
- } else {
+ } else {
mindlinContactPhysics->Delta1 = max(solution.delta1,solution.delta2);
- mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2);}
+ mindlinContactPhysics->Delta2 = min(solution.delta1,solution.delta2);
+ }
}
- else if (fusionDetection) bodiesMenisciiList.remove((*ii));//If the interaction is not real, it should not be in the list
- }
- }
- }
-
- if (fusionDetection) checkFusion();
-
- for(ii= scene->interactions->begin(); ii!=iiEnd ; ++ii )
+// }
+ } else /*not real*/ if (fusionDetection) bodiesMenisciiList.remove((*ii));//If the interaction is not real, it should not be in the list
+ }
+ if (fusionDetection) checkFusion();
+
+ for (ii= scene->interactions->begin(); ii!=iiEnd ; ++ii)
{
if ((*ii)->isReal())
{
@@ -193,15 +182,15 @@
MindlinCapillaryPhys* mindlinContactPhysics=NULL;
if (!hertzOn) cundallContactPhysics = static_cast<CapillaryPhys*>((*ii)->phys.get());//use CapillaryPhys for linear model
else mindlinContactPhysics = static_cast<MindlinCapillaryPhys*>((*ii)->phys.get());//use MindlinCapillaryPhys for hertz model
-
- if ((hertzOn && mindlinContactPhysics->meniscus) || (!hertzOn && cundallContactPhysics->meniscus))
+
+ if ((hertzOn && mindlinContactPhysics->meniscus) || (!hertzOn && cundallContactPhysics->meniscus))
{
if (fusionDetection)
{//version with effect of fusion
//BINARY VERSION : if fusionNumber!=0 then no capillary force
short int& fusionNumber = hertzOn?mindlinContactPhysics->fusionNumber:cundallContactPhysics->fusionNumber;
if (binaryFusion)
- {
+ {
if (fusionNumber!=0)
{ //cerr << "fusion" << endl;
hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap = Vector3r::Zero();
@@ -216,7 +205,7 @@
scene->forces.addForce((*ii)->getId2(),-(hertzOn?mindlinContactPhysics->Fcap:cundallContactPhysics->Fcap));
}
}
- }
+ }
}
capillarylaw::capillarylaw()
@@ -243,10 +232,8 @@
if ( !bodiesMenisciiList[i].empty() )
{
lastMeniscus = bodiesMenisciiList[i].end();
- //cerr << "size = "<<bodiesMenisciiList[i].size() << " empty="<<bodiesMenisciiList[i].empty() <<endl;
for ( firstMeniscus=bodiesMenisciiList[i].begin(); firstMeniscus!=lastMeniscus; ++firstMeniscus )//FOR EACH MENISCUS ON THIS BODY...
{
- //if (*firstMeniscus)->isReal();
CapillaryPhys* interactionPhysics1 = YADE_CAST<CapillaryPhys*>((*firstMeniscus)->phys.get());
currentMeniscus = firstMeniscus; ++currentMeniscus;
@@ -266,23 +253,15 @@
Vector3r normalFirstMeniscus = YADE_CAST<ScGeom*>((*firstMeniscus)->geom.get())->normal;
Vector3r normalCurrentMeniscus = YADE_CAST<ScGeom*>((*currentMeniscus)->geom.get())->normal;
- //if (i != (*firstMeniscus)->getId1()) normalFirstMeniscus = -normalFirstMeniscus;
- //if (i != (*currentMeniscus)->getId1()) normalCurrentMeniscus = -normalCurrentMeniscus;
- // normal is always from id1 to id2
-
Real normalDot = 0;
if ((*firstMeniscus)->getId1() == (*currentMeniscus)->getId1() || (*firstMeniscus)->getId2() == (*currentMeniscus)->getId2()) normalDot = normalFirstMeniscus.dot(normalCurrentMeniscus);
else
normalDot = - (normalFirstMeniscus.dot(normalCurrentMeniscus));
- //cerr << "normalDot ="<< normalDot << endl;
-
Real normalAngle = 0;
if (normalDot >= 0 ) normalAngle = Mathr::FastInvCos0(normalDot);
else normalAngle = ((Mathr::PI) - Mathr::FastInvCos0(-(normalDot)));
- //cerr << "sumMeniscAngle= " << (angle1+angle2)<< "| normalAngle" << normalAngle*Mathr::RAD_TO_DEG << endl;
-
if ((angle1+angle2)*Mathr::DEG_TO_RAD > normalAngle)
{
++(interactionPhysics1->fusionNumber); ++(interactionPhysics2->fusionNumber);//count +1 if 2 meniscii are overlaping
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2011-10-06 15:31:48 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2011-11-30 17:39:33 +0000
@@ -105,7 +105,7 @@
#endif
YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
// clump members are handled inside clumps
- if(unlikely(!b || b->isClumpMember())) continue;
+ if(unlikely(b->isClumpMember())) continue;
State* state=b->state.get(); const Body::id_t& id=b->getId();
Vector3r f=gravity*state->mass, m=Vector3r::Zero();
