------------------------------------------------------------
revno: 3590
committer: Christian Jakob <ja...@ifgt.tu-freiberg.de>
timestamp: Thu 2015-02-19 13:54:58 +0100
message:
fix fusion detection for hertz model in capillary law
modified:
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to
https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2014-12-10 01:37:09 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2015-02-19 12:54:58 +0000
@@ -58,18 +58,25 @@
if (!scene) cerr << "scene not defined!";
if (!capillary) postLoad(*this);//when the script does not define arguments, postLoad is never called
shared_ptr<BodyContainer>& bodies = scene->bodies;
- if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene);
+
+ //check for contact model once (assuming that contact model does not change)
+ if (!hertzInitialized){//NOTE: We are assuming that only one type is used in one simulation here
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if (I->isReal()) {
+ if (CapillaryPhys::getClassIndexStatic()==I->phys->getClassIndex()) hertzOn=false;
+ else if (MindlinCapillaryPhys::getClassIndexStatic()==I->phys->getClassIndex()) hertzOn=true;
+ else LOG_ERROR("The capillary law is not implemented for interactions using"<<I->phys->getClassName());
+ break;
+ }
+ }
+ hertzInitialized = true;
+ }
+
+ if (fusionDetection && !bodiesMenisciiList.initialized) bodiesMenisciiList.prepare(scene,hertzOn);
- bool hertzInitialized = false;
FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){ // could be done in parallel ? Was tried once, but needs maybe more comparison to assert it is OK: http://www.mail-archive.com/yade-dev@lists.launchpad.net/msg10842.html
/// interaction is real
if (interaction->isReal()) {
- 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());
- }
- hertzInitialized = true;
CapillaryPhys* cundallContactPhysics=NULL;
MindlinCapillaryPhys* mindlinContactPhysics=NULL;
@@ -148,7 +155,7 @@
if (!Vinterpol) {
if ((fusionDetection) || (hertzOn ? mindlinContactPhysics->isBroken : cundallContactPhysics->isBroken)) bodiesMenisciiList.remove(interaction);
if (D>0) scene->interactions->requestErase(interaction);
- else if (Pinterpol > 0) LOG_ERROR("No meniscus found at a contact. capillaryPressure may be too large wrt. the loaded data files.") // V=0 at a contact reveals a problem if and only if uc* > 0
+ else if ((Pinterpol > 0)) LOG_ERROR("No meniscus found at a contact. capillaryPressure may be too large wrt. the loaded data files.") // V=0 at a contact reveals a problem if and only if uc* > 0
}
/// wetting angles
if (!hertzOn) {
@@ -204,9 +211,10 @@
{
//Reset fusion numbers
FOREACH(const shared_ptr<Interaction>& interaction, *scene->interactions){ // same remark for parallel loops, the most problematic part ?
- if ( !interaction->isReal()) continue;
- if (!hertzOn) static_cast<CapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
- else static_cast<MindlinCapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
+ if ( interaction->isReal()) {
+ if (!hertzOn) static_cast<CapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
+ else static_cast<MindlinCapillaryPhys*>(interaction->phys.get())->fusionNumber=0;
+ }
}
list< shared_ptr<Interaction> >::iterator firstMeniscus, lastMeniscus, currentMeniscus;
@@ -260,7 +268,7 @@
if ((angle1+angle2)*Mathr::DEG_TO_RAD > normalAngle) {
if (!hertzOn) {++(cundallInteractionPhysics1->fusionNumber); ++(cundallInteractionPhysics2->fusionNumber);}//count +1 if 2 meniscii are overlaping
else {++(mindlinInteractionPhysics1->fusionNumber); ++(mindlinInteractionPhysics2->fusionNumber);}
- };
+ }
}
}
}
@@ -499,13 +507,13 @@
return os;
}
-BodiesMenisciiList::BodiesMenisciiList(Scene * scene)
+BodiesMenisciiList::BodiesMenisciiList(Scene * scene, bool hertzOn)
{
initialized=false;
- prepare(scene);
+ prepare(scene, hertzOn);
}
-bool BodiesMenisciiList::prepare(Scene * scene)
+bool BodiesMenisciiList::prepare(Scene * scene, bool hertzOn)
{
//cerr << "preparing bodiesInteractionsList" << endl;
interactionsOnBody.clear();
@@ -523,10 +531,11 @@
{
interactionsOnBody[i].clear();
}
-
+
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){ // parallel version using Engine::ompThreads variable not accessible, this function is not one of the Engine..
if (I->isReal()) {
- if (static_cast<CapillaryPhys*>(I->phys.get())->meniscus) insert(I);
+ if (!hertzOn) {if (static_cast<CapillaryPhys*>(I->phys.get())->meniscus) insert(I);}
+ else {if (static_cast<MindlinCapillaryPhys*>(I->phys.get())->meniscus) insert(I);}
}
}
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2015-02-03 21:52:41 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2015-02-19 12:54:58 +0000
@@ -48,7 +48,7 @@
/// R = ratio(RadiusParticle1 on RadiusParticle2). Here, 10 R values from interpolation files (yade/extra/capillaryFiles), R = 1, 1.1, 1.25, 1.5, 1.75, 2, 3, 4, 5, 10
const int NB_R_VALUES = 10;
-class capillarylaw; // fait appel a la classe def plus bas
+class capillarylaw; // fait appel a la classe def plus bas //TODO: translate this in english
class Interaction;
///This container class is used to check if meniscii overlap. Wet interactions are put in a series of lists, with one list per body.
@@ -61,8 +61,8 @@
public:
BodiesMenisciiList();
- BodiesMenisciiList(Scene*);
- bool prepare(Scene*);
+ BodiesMenisciiList(Scene*,bool);
+ bool prepare(Scene*,bool);
bool insert(const shared_ptr<Interaction>&);
bool remove(const shared_ptr<Interaction>&);
list< shared_ptr<Interaction> >& operator[] (int);
@@ -70,7 +70,6 @@
void display();
void checkLengthBuffer(const shared_ptr<Interaction>&);
-
bool initialized;
};
@@ -85,15 +84,21 @@
void action();
void postLoad(Law2_ScGeom_CapillaryPhys_Capillarity&);
+ bool hertzInitialized;
+ bool hertzOn;
+ bool getHertzOn();
+
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(Law2_ScGeom_CapillaryPhys_Capillarity,GlobalEngine,"This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).\n\nThe control parameter is the :yref:`capillary pressure<Law2_ScGeom_CapillaryPhys_Capillarity::capillaryPressure>` (or suction) Uc = Ugas - Uliquid. Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of the defined capillary pressure and of the interacting geometry (spheres radii and interparticular distance).\n\nReferences: in english [Scholtes2009b]_; more detailed, but in french [Scholtes2009d]_.\n\nThe law needs ascii files M(r=i) with i=R1/R2 to work (see https://yade-dem.org/wiki/CapillaryTriaxialTest). These ASCII files contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry, assuming a null wetting angle.\n\nIn order to allow capillary forces between distant spheres, it is necessary to enlarge the bounding boxes using :yref:`Bo1_Sphere_Aabb::aabbEnlargeFactor` and make the Ig2 define define distant interactions via :yref:`interactionDetectionFactor<Ig2_Sphere_Sphere_ScGeom::interactionDetectionFactor>`. It is also necessary to disable interactions removal by the constitutive law (:yref:`Law2<Law2_ScGeom_FrictPhys_CundallStrack::neverErase>=True`). The only combinations of laws supported are currently capillary law + :yref:`Law2_ScGeom_FrictPhys_CundallStrack` and capillary law + :yref:`Law2_ScGeom_MindlinPhys_Mindlin` (and the other variants of Hertz-Mindlin).\n\nSee CapillaryPhys-example.py for an example script.",
((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defined as Uc=Ugas-Uliquid"))
((bool,fusionDetection,false,,"If true potential menisci overlaps are checked, computing :yref:`fusionNumber<CapillaryPhys.fusionNumber>` for each capillary interaction, and reducing :yref:`fCap<CapillaryPhys.fCap>` according to :yref:`binaryFusion<Law2_ScGeom_CapillaryPhys_Capillarity.binaryFusion>`"))
((bool,binaryFusion,true,,"If true, capillary forces are set to zero as soon as, at least, 1 overlap (menisci fusion) is detected. Otherwise :yref:`fCap<CapillaryPhys.fCap>` = :yref:`fCap<CapillaryPhys.fCap>` / (:yref:`fusionNumber<CapillaryPhys.fusionNumber>` + 1 )"))
- ((bool,hertzOn,false,,"|yupdate| true if hertz model is used"))
((bool,createDistantMeniscii,false,,"Generate meniscii between distant spheres? Else only maintain the existing ones. For modeling a wetting path this flag should always be false. For a drying path it should be true for one step (initialization) then false, as in the logic of [Scholtes2009c]_"))
,/*deprec*/
((CapillaryPressure,capillaryPressure,"naming convention"))
- ,,,
+ ,,/*constructor*/
+ hertzInitialized = false;
+ hertzOn = false;
+ ,
);
};
@@ -108,7 +113,7 @@
~TableauD();
};
-// Fonction d'ecriture de tableau, utilisee dans le constructeur pour test
+// Fonction d'ecriture de tableau, utilisee dans le constructeur pour test //TODO: translate this in english
class Tableau;
std::ostream& operator<<(std::ostream& os, Tableau& T);
_______________________________________________
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
