yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10981
[Branch ~yade-pkg/yade/git-trunk] Rev 4030: Liquid migration model, code refactoring.
------------------------------------------------------------
revno: 4030
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-06-18 09:35:57 +0200
message:
Liquid migration model, code refactoring.
modified:
examples/capillary/liquidmigration/showcase.py
pkg/dem/ViscoelasticCapillarPM.cpp
--
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 'examples/capillary/liquidmigration/showcase.py'
--- examples/capillary/liquidmigration/showcase.py 2014-05-08 05:57:04 +0000
+++ examples/capillary/liquidmigration/showcase.py 2014-06-18 07:35:57 +0000
@@ -29,8 +29,8 @@
id5 = O.bodies.append(sphere(center=[(r1+r2)*d*2,(r1+r2)*d,0], radius=r2,material=mat1,fixed=True, color=[0,1,0]))
-Vf = 0.5e-1
-Vfmin = 0.1e-1
+Vf = 0.0e-1
+Vfmin = 0.0e-1
O.bodies[id1].Vf = Vf
O.bodies[id1].Vmin = Vfmin
@@ -61,7 +61,7 @@
[Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],
[Law2_ScGeom_ViscElCapPhys_Basic()],
),
- LiqControl(),
+ LiqControl(label='lqc'),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
PyRunner(command='showData()',iterPeriod=1),
]
@@ -105,14 +105,18 @@
O.bodies[i].Vf = 0
O.bodies[i].Vmin = 0
-O.interactions[id1,id2].phys.Vb = 1.0
O.interactions[id1,id2].phys.Vmax = 5.0
-O.interactions[id2,id3].phys.Vb = 1.0
+lqc.addLiqInter(id1, id2, 1.0)
+
O.interactions[id2,id3].phys.Vmax = 5.0
-O.interactions[id3,id4].phys.Vb = 1.0
+lqc.addLiqInter(id2, id3, 1.0)
+
O.interactions[id3,id4].phys.Vmax = 5.0
-O.interactions[id3,id5].phys.Vb = 1.0
+lqc.addLiqInter(id3, id4, 1.0)
+
O.interactions[id3,id5].phys.Vmax = 5.0
+lqc.addLiqInter(id3, id5, 1.0)
+
O.run(1, True)
=== modified file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 2014-06-13 13:58:12 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-06-18 07:35:57 +0000
@@ -361,8 +361,13 @@
// Calculate, how much new contacts will be at each body
for (unsigned int i=0; i<scene->addIntrs.size(); i++) {
- addBodyMapInt( bI, scene->addIntrs[i]->getId1() );
- addBodyMapInt( bI, scene->addIntrs[i]->getId2() );
+ shared_ptr<Body> b1 = Body::byId(scene->addIntrs[i]->getId1(),scene);
+ shared_ptr<Body> b2 = Body::byId(scene->addIntrs[i]->getId2(),scene);
+
+ if(not(mask!=0 && ((b1->groupMask & b2->groupMask & mask)==0))) {
+ addBodyMapInt( bI, scene->addIntrs[i]->getId1() );
+ addBodyMapInt( bI, scene->addIntrs[i]->getId2() );
+ }
}
// Update volume water at each deleted interaction for each body
@@ -384,43 +389,42 @@
const id_t id2 = b2->id;
ViscElCapPhys* Vb=dynamic_cast<ViscElCapPhys*>(scene->addIntrs[i]->phys.get());
- const Real Vmax = vMax(b1, b2);
- Vb->Vmax = Vmax;
-
- Real Vf1 = 0.0;
- Real Vf2 = 0.0;
-
- if ((b1->Vmin)<b1->Vf) {
- Vf1 = (b1->Vf - b1->Vmin)/bI[id1];
- }
-
- if ((b2->Vmin)<b2->Vf) {
- Vf2 = (b2->Vf - b2->Vmin)/bI[id2];
- }
-
- Real Vrup = Vf1+Vf2;
-
+
if(mask!=0 && ((b1->groupMask & b2->groupMask & mask)==0)) {
- Vf1 = 0;
- Vf2 = 0;
- Vrup = 0;
- } else if (Vrup > Vmax) {
- Vf1 *= Vmax/Vrup;
- Vf2 *= Vmax/Vrup;
- Vrup = Vf1 + Vf2;
- }
-
- liqVolShr += Vrup;
- addBodyMapReal(bodyUpdateLiquid, id1, -Vf1);
- addBodyMapReal(bodyUpdateLiquid, id2, -Vf2);
-
- Vb->Vb = Vrup;
- if (particleconserve) {
- Vb->Vf1 = Vf1;
- Vb->Vf2 = Vf2;
+ Vb->Vb = 0.0;
+ Vb->Vf1 = 0.0;
+ Vb->Vf2 = 0.0;
+ Vb->Vmax = 0.0;
} else {
- Vb->Vf1 = Vrup/2.0;
- Vb->Vf2 = Vrup/2.0;
+ const Real Vmax = vMax(b1, b2);
+ Vb->Vmax = Vmax;
+
+ Real Vf1 = 0.0;
+ Real Vf2 = 0.0;
+
+ if ((b1->Vmin)<b1->Vf) { Vf1 = (b1->Vf - b1->Vmin)/bI[id1]; }
+ if ((b2->Vmin)<b2->Vf) { Vf2 = (b2->Vf - b2->Vmin)/bI[id2]; }
+
+ Real Vrup = Vf1+Vf2;
+
+ if (Vrup > Vmax) {
+ Vf1 *= Vmax/Vrup;
+ Vf2 *= Vmax/Vrup;
+ Vrup = Vf1 + Vf2;
+ }
+
+ liqVolShr += Vrup;
+ addBodyMapReal(bodyUpdateLiquid, id1, -Vf1);
+ addBodyMapReal(bodyUpdateLiquid, id2, -Vf2);
+
+ Vb->Vb = Vrup;
+ if (particleconserve) {
+ Vb->Vf1 = Vf1;
+ Vb->Vf2 = Vf2;
+ } else {
+ Vb->Vf1 = Vrup/2.0;
+ Vb->Vf2 = Vrup/2.0;
+ }
}
}
@@ -435,7 +439,6 @@
for (mapBodyInt::const_iterator it = bodyNeedUpdate.begin(); it != bodyNeedUpdate.end(); ++it) {
updateLiquid(Body::byId(it->first));
}
-
}
void LiqControl::updateLiquid(shared_ptr<Body> b){
@@ -451,7 +454,7 @@
for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
if(!((*it).second) or !(((*it).second)->isReal())) continue;
ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
- if (physT->Vb<physT->Vmax) {
+ if ((physT->Vb < physT->Vmax) and (physT->Vmax > 0)) {
LiqContactsAccept+=physT->Vmax-physT->Vb;
contactN++;
}
@@ -472,15 +475,15 @@
for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
if(!((*it).second) or !(((*it).second)->isReal())) continue;
ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
- if (physT->Vb<physT->Vmax) {
+ if ((physT->Vb < physT->Vmax) and (physT->Vmax > 0)) {
const Real addVolLiq = (physT->Vmax - physT->Vb)*FillLevel;
liqVolShr += addVolLiq;
physT->Vb += addVolLiq;
if (particleconserve) {
if (((*it).second)->getId1() == (*it).first) {
+ physT->Vf2+=addVolLiq;
+ } else if (((*it).second)->getId2() == (*it).first) {
physT->Vf1+=addVolLiq;
- } else if (((*it).second)->getId2() == (*it).first) {
- physT->Vf2+=addVolLiq;
}
} else {
physT->Vf1+=addVolLiq/2.0;
@@ -538,7 +541,19 @@
if(!((*it).second) or !(((*it).second)->isReal())) continue;
ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
if (physT and physT->Vb and physT->Vb>0) {
- LiqVol += physT->Vb/2.0;
+ if (((*it).second)->id1 == b->id) {
+ if (physT->Vf1 > 0 or physT->Vf2 > 0) {
+ LiqVol += physT->Vf1;
+ } else {
+ LiqVol += physT->Vb/2.0;
+ }
+ } else {
+ if (physT->Vf1 > 0 or physT->Vf2 > 0) {
+ LiqVol += physT->Vf2;
+ } else {
+ LiqVol += physT->Vb/2.0;
+ }
+ }
}
}
return LiqVol;
