yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11456
[Branch ~yade-pkg/yade/git-trunk] Rev 3452: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Christian Jakob (jakob-ifgt)
Jan Stránský (honzik)...
------------------------------------------------------------
revno: 3452 [merge]
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Wed 2014-10-08 15:01:56 +0200
message:
Merge github.com:yade/trunk into chaoUnsat
removed:
pkg/pfv/FlowEngine.ipp
modified:
core/InteractionContainer.hpp
doc/references.bib
doc/sphinx/github.rst
doc/sphinx/introduction.rst
doc/sphinx/prog.rst
doc/sphinx/yadeSphinx.py
doc/yade-articles.bib
examples/polyhedra/ball.py
examples/polyhedra/free-fall.py
examples/polyhedra/irregular.py
examples/polyhedra/wall.py
examples/tetra/twoTetrasPoly.py
examples/triax-tutorial/script-session1.py
gui/qt4/SerializableEditor.py
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/FlowBoundingSphereLinSolv.ipp
lib/triangulation/PeriodicFlow.hpp
pkg/common/GravityEngines.cpp
pkg/common/Grid.cpp
pkg/common/InsertionSortCollider.cpp
pkg/common/PeriodicEngines.hpp
pkg/dem/CapillaryPhys.hpp
pkg/dem/CohesiveFrictionalContactLaw.hpp
pkg/dem/HertzMindlin.hpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp
pkg/dem/NewtonIntegrator.cpp
pkg/dem/Polyhedra.cpp
pkg/dem/Polyhedra.hpp
pkg/dem/Polyhedra_Ig2.cpp
pkg/dem/Polyhedra_Ig2.hpp
pkg/dem/SpherePack.cpp
pkg/dem/Tetra.cpp
pkg/dem/TriaxialStressController.cpp
pkg/dem/VTKRecorder.cpp
pkg/dem/ViscoelasticCapillarPM.cpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.hpp.in
pkg/pfv/FlowEngine.ipp.in
py/post2d.py
py/utils.py
py/wrapper/yadeWrapper.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 'core/InteractionContainer.hpp'
--- core/InteractionContainer.hpp 2014-07-03 17:20:40 +0000
+++ core/InteractionContainer.hpp 2014-09-15 10:27:11 +0000
@@ -103,12 +103,35 @@
*/
template<class T> size_t conditionalyEraseNonReal(const T& t, Scene* rb){
// beware iterators here, since erase is invalidating them. We need to iterate carefully, and keep in mind that erasing one interaction is moving the last one to the current position.
- size_t initSize=currSize;
- for (size_t linPos=0; linPos<currSize;){
- const shared_ptr<Interaction>& i=linIntrs[linPos];
- if(!i->isReal() && t.shouldBeErased(i->getId1(),i->getId2(),rb)) erase(i->getId1(),i->getId2(),linPos);
- else linPos++;}
- return initSize-currSize;
+ // For the parallel flavor we build the list to be erased in parallel, then it is erased sequentially. Still significant speedup since checking bounds is the most expensive part.
+ #ifdef YADE_OPENMP
+ if (omp_get_max_threads()<=1) {
+ #endif
+ size_t initSize=currSize;
+ for (size_t linPos=0; linPos<currSize;){
+ const shared_ptr<Interaction>& i=linIntrs[linPos];
+ if(!i->isReal() && t.shouldBeErased(i->getId1(),i->getId2(),rb)) erase(i->getId1(),i->getId2(),linPos);
+ else linPos++;}
+ return initSize-currSize;
+ #ifdef YADE_OPENMP
+ } else {
+ unsigned nThreads= omp_get_max_threads();
+ assert(nThreads>0);
+ std::vector<std::vector<Vector3i > >toErase;
+ toErase.resize(nThreads,std::vector<Vector3i >());
+ for (unsigned kk=0; kk<nThreads; kk++) toErase[kk].reserve(1000);//A smarter value than 1000?
+ size_t initSize=currSize;
+ #pragma omp parallel for schedule(static) num_threads(nThreads)
+ for (size_t linPos=0; linPos<currSize;linPos++){
+ const shared_ptr<Interaction>& i=linIntrs[linPos];
+ if(!i->isReal() && t.shouldBeErased(i->getId1(),i->getId2(),rb)) toErase[omp_get_thread_num()].push_back(Vector3i(i->getId1(),i->getId2(),linPos)) ;
+ }
+ for (int kk=nThreads-1; kk>=0; kk--)
+ for (int ii(toErase[kk].size()-1); ii>=0;ii--)
+ erase(toErase[kk][ii][0],toErase[kk][ii][1],toErase[kk][ii][2]);
+ return initSize-currSize;
+ }
+ #endif
}
// we must call Scene's ctor (and from Scene::postLoad), since we depend on the existing BodyContainer at that point.
=== modified file 'doc/references.bib'
--- doc/references.bib 2014-08-19 13:59:12 +0000
+++ doc/references.bib 2014-09-12 05:57:43 +0000
@@ -670,7 +670,7 @@
number={6},
doi={10.1007/s10035-007-0065-z},
title={Coefficient of restitution and linear–dashpot model revisited},
- url={http://dx.doi.org/10.1007/s10035-007-0065-z},
+ url={http://arxiv.org/pdf/cond-mat/0701278},
publisher={Springer-Verlag},
keywords={Particle collisions; Coefficient of restitution},
author={Schwager, Thomas and Pöschel, Thorsten},
=== modified file 'doc/sphinx/github.rst'
--- doc/sphinx/github.rst 2014-08-18 15:34:10 +0000
+++ doc/sphinx/github.rst 2014-10-08 05:24:38 +0000
@@ -150,7 +150,9 @@
[branch "master"]
remote = origin
+
merge = refs/heads/master
+
rebase = true
Auto-rebase may have unpleasant side effects by blocking "pull" if you have uncommited changes. In this case you can use "git stash"::
=== modified file 'doc/sphinx/introduction.rst'
--- doc/sphinx/introduction.rst 2013-09-27 12:13:27 +0000
+++ doc/sphinx/introduction.rst 2014-10-08 05:24:38 +0000
@@ -131,7 +131,7 @@
Running simulation
------------------
-As explained above, the loop consists in running defined sequence of engines. Step number can be queried by ``O.iter`` and advancing by one step is done by ``O.step()``. Every step advances *virtual time* by current timestep, ``O.dt``:
+As explained below, the loop consists in running defined sequence of engines. Step number can be queried by ``O.iter`` and advancing by one step is done by ``O.step()``. Every step advances *virtual time* by current timestep, ``O.dt``:
.. ipython::
=== modified file 'doc/sphinx/prog.rst'
--- doc/sphinx/prog.rst 2014-08-21 08:51:47 +0000
+++ doc/sphinx/prog.rst 2014-09-18 09:28:42 +0000
@@ -1508,14 +1508,10 @@
#. For real interactions, :yref:`InteractionLoop` (via :yref:`LawDispatcher`) calls appropriate :yref:`LawFunctor` based on combination of :yref:`IGeom` and :yref:`IPhys` of the interaction. Again, it is an error if no functor capable of handling it is found.
-#. :yref:`LawDispatcher` can decide that an interaction should be removed (such as if bodies get too far apart for non-cohesive laws; or in case of complete damage for damage models). This is done by calling
-
- .. code-block:: c++
-
- InteractionContainer::requestErase(id1,id2)
-
- Such interaction will not be deleted immediately, but will be reset to potential state. At next step, the collider will call ``InteractionContainer::erasePending``, which will only completely erase interactions the collider indicates; the rest will be kept in potential state.
-
+#. :yref:`LawDispatcher` takes care of erasing those interactions that are no longer active (such as if bodies get too far apart for non-cohesive laws; or in case of complete damage for damage models). This is triggered by the :yref:`LawFunctor` returning false. For this reason it is of upmost importance for the :yref:`LawFunctor` to return consistently.
+
+Such interaction will not be deleted immediately, but will be reset to potential state. At the next execution of the collider ``InteractionContainer::conditionalyEraseNonReal`` will be called, which will completely erase interactions only if the bounding boxes ceased to overlap; the rest will be kept in potential state.
+
Creating interactions explicitly
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
=== modified file 'doc/sphinx/yadeSphinx.py'
--- doc/sphinx/yadeSphinx.py 2013-10-15 16:14:46 +0000
+++ doc/sphinx/yadeSphinx.py 2014-10-02 08:56:07 +0000
@@ -212,12 +212,10 @@
global writer
writer=None
-for w in ['html','latex','epub']:
- writer=w
+for writer in ['html','latex','epub']:
genWrapperRst()
- # HACK: must rewrite sys.argv, since reference generator in conf.py determines if we output latex/html by inspecting it
- sys.argv=['sphinx-build','-a','-E','-b','%s'%writer,'-d',outDir+'/doctrees','.',outDir+'/%s'%writer]
- sphinx.main(sys.argv)
+ commLine = "sphinx-build -a -E -b %s -d %s . %s"%(writer,outDir+'/doctrees', (outDir+'/%s'%writer))
+ os.system(commLine)
if writer=='html':
makeBaseClassesClickable((outDir+'/html/yade.wrapper.html'),writer)
elif writer=='latex':
@@ -225,7 +223,7 @@
if (os.path.exists('/usr/share/javascript/jquery/jquery.js')): #Check, whether jquery.js installed in system
os.system('rm '+ outDir+'/html/_static/jquery.js')
os.system('ln -s /usr/share/javascript/jquery/jquery.js '+ outDir+'/html/_static/jquery.js')
-
+
# HACK!!!!==========================================================================
# New sphinx-python versions (hopefully) are producing empty "verbatim"-environments.
# That is why xelatex crashes.
=== modified file 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2014-04-02 15:19:41 +0000
+++ doc/yade-articles.bib 2014-09-12 15:26:09 +0000
@@ -604,3 +604,17 @@
doi = {10.1103/PhysRevE.88.062203},
url = {http://link.aps.org/doi/10.1103/PhysRevE.88.062203}
}
+
+@article{Gladky2014,
+ year={2014},
+ issn={1434-5021},
+ journal={Granular Matter},
+ doi={10.1007/s10035-014-0527-z},
+ title={Comparison of different capillary bridge models for application in the discrete element method},
+ url={http://dx.doi.org/10.1007/s10035-014-0527-z},
+ publisher={Springer Berlin Heidelberg},
+ keywords={Granular material; DEM; Capillary bridge models; Liquid bridges},
+ author={Gladkyy, Anton and Schwarze, Rüdiger},
+ pages={1-10},
+ language={English}
+}
=== modified file 'examples/polyhedra/ball.py'
--- examples/polyhedra/ball.py 2013-10-16 15:24:49 +0000
+++ examples/polyhedra/ball.py 2014-09-22 12:28:50 +0000
@@ -27,7 +27,7 @@
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
- [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] # contact law -- apply forces
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(damping=0.5,gravity=(0,0,-9.81)),
=== modified file 'examples/polyhedra/free-fall.py'
--- examples/polyhedra/free-fall.py 2013-10-16 15:24:49 +0000
+++ examples/polyhedra/free-fall.py 2014-09-22 12:28:50 +0000
@@ -36,7 +36,7 @@
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
- [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] # contact law -- apply forces
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(damping=0.3,gravity=(0,0,-9.81)),
=== modified file 'examples/polyhedra/irregular.py'
--- examples/polyhedra/irregular.py 2013-10-16 15:24:49 +0000
+++ examples/polyhedra/irregular.py 2014-09-22 12:28:50 +0000
@@ -38,7 +38,7 @@
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
- [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] # contact law -- apply forces
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(damping=0.5,gravity=(0,0,-9.81)),
=== modified file 'examples/polyhedra/wall.py'
--- examples/polyhedra/wall.py 2013-10-16 15:24:49 +0000
+++ examples/polyhedra/wall.py 2014-09-22 12:28:50 +0000
@@ -50,7 +50,7 @@
InteractionLoop(
[Ig2_Wall_Polyhedra_PolyhedraGeom(), Ig2_Polyhedra_Polyhedra_PolyhedraGeom(), Ig2_Facet_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()], # collision "physics"
- [PolyhedraVolumetricLaw()] # contact law -- apply forces
+ [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()] # contact law -- apply forces
),
#GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(damping=0.5,gravity=(0,0,-9.81)),
=== modified file 'examples/tetra/twoTetrasPoly.py'
--- examples/tetra/twoTetrasPoly.py 2013-10-15 17:11:37 +0000
+++ examples/tetra/twoTetrasPoly.py 2014-09-22 12:28:50 +0000
@@ -27,7 +27,7 @@
InteractionLoop(
[Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
[Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys()],
- [PolyhedraVolumetricLaw()]
+ [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
),
NewtonIntegrator(damping=0),
PyRunner(iterPeriod=500,command="addPlotData()"),
=== modified file 'examples/triax-tutorial/script-session1.py'
--- examples/triax-tutorial/script-session1.py 2014-08-01 15:37:51 +0000
+++ examples/triax-tutorial/script-session1.py 2014-09-22 23:13:24 +0000
@@ -80,7 +80,7 @@
############################
triax=TriaxialStressController(
- ## ThreeDTriaxialEngine will be used to control stress and strain. It controls particles size and plates positions.
+ ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
=== modified file 'gui/qt4/SerializableEditor.py'
--- gui/qt4/SerializableEditor.py 2014-06-12 16:54:31 +0000
+++ gui/qt4/SerializableEditor.py 2014-10-01 17:32:47 +0000
@@ -18,14 +18,11 @@
seqSerializableShowType=True # show type headings in serializable sequences (takes vertical space, but makes the type hyperlinked)
-# BUG: cursor is moved to the beginnign of the input field even if it has focus
-#
+# BUG: cursor is moved to the beginning of the input field even if it has focus
# checking for focus seems to return True always and cursor is never moved
-#
# the 'True or' part effectively disables the condition (so that the cursor is moved always), but it might be fixed in the future somehow
-#
# if True or w.hasFocus(): w.home(False)
-#
+# Janek: It looks like I've fixed this BUG, please do more testing.
#
def makeWrapperHref(text,className,attr=None,static=False):
@@ -101,7 +98,8 @@
QSpinBox.__init__(self,parent)
self.setRange(int(-1e9),int(1e9)); self.setSingleStep(1);
self.valueChanged.connect(self.update)
- def refresh(self): self.setValue(self.getter())
+ def refresh(self):
+ if (not self.hasFocus()): self.setValue(self.getter())
def update(self): self.trySetter(self.value())
class AttrEditor_Str(AttrEditor,QLineEdit):
@@ -111,7 +109,8 @@
self.textEdited.connect(self.isHot)
self.selectionChanged.connect(self.isHot)
self.editingFinished.connect(self.update)
- def refresh(self): self.setText(self.getter())
+ def refresh(self):
+ if (not self.hasFocus()): self.setText(self.getter())
def update(self): self.trySetter(str(self.text()))
class AttrEditor_Float(AttrEditor,QLineEdit):
@@ -122,12 +121,53 @@
self.selectionChanged.connect(self.isHot)
self.editingFinished.connect(self.update)
def refresh(self):
- self.setText(str(self.getter()));
- if True or not self.hasFocus(): self.home(False)
+ #if True or not self.hasFocus(): self.home(False)
+ if (not self.hasFocus()):
+ self.setText(str(self.getter()));
+ self.home(False)
def update(self):
try: self.trySetter(float(self.text()))
except ValueError: self.refresh()
+class AttrEditor_Complex(AttrEditor,QLineEdit):
+ def __init__(self,parent,getter,setter):
+ AttrEditor.__init__(self,getter,setter)
+ QFrame.__init__(self,parent)
+ self.rows,self.cols=1,2
+ self.setContentsMargins(0,0,0,0)
+ self.first=True
+ val=self.getter()
+ self.grid=QGridLayout(self); self.grid.setSpacing(0); self.grid.setMargin(0)
+ for row,col in itertools.product(range(self.rows),range(self.cols)):
+ w=QLineEdit('')
+ self.grid.addWidget(w,row,col);
+ w.textEdited.connect(self.isHot)
+ w.selectionChanged.connect(self.isHot)
+ w.editingFinished.connect(self.update)
+ def refresh(self,force=False):
+ val=self.getter()
+ for row,col in itertools.product(range(self.rows),range(self.cols)):
+ w=self.grid.itemAtPosition(row,col).widget()
+ if(self.first or force):
+ w.setText(str(val.real if col==0 else val.imag))
+ #if True or not w.hasFocus: w.home(False) # make the left-most part visible, if the text is wider than the widget
+ if (not w.hasFocus):
+ w.setText(str(val.real if col==0 else val.imag))
+ w.home(False) # make the left-most part visible, if the text is wider than the widget
+ self.first=False
+ def update(self):
+ try:
+ val=self.getter()
+ w1=self.grid.itemAtPosition(0,0).widget()
+ w2=self.grid.itemAtPosition(0,1).widget()
+ if w1.isModified() or w2.isModified():
+ val=complex(float(w1.text()),float(w2.text()))
+ logging.debug('setting'+str(val))
+ self.trySetter(val)
+ except ValueError:
+ self.refresh(force=True)
+ def setFocus(self): self.grid.itemAtPosition(0,0).widget().setFocus()
+
class AttrEditor_Quaternion(AttrEditor,QFrame):
def __init__(self,parent,getter,setter):
AttrEditor.__init__(self,getter,setter)
@@ -135,7 +175,8 @@
self.grid=QHBoxLayout(self); self.grid.setSpacing(0); self.grid.setMargin(0)
for i in range(4):
if i==3:
- f=QFrame(self); f.setFrameShape(QFrame.VLine); f.setFrameShadow(QFrame.Sunken); f.setFixedWidth(4) # add vertical divider (axis | angle)
+ # add vertical divider (axis | angle)
+ f=QFrame(self); f.setFrameShape(QFrame.VLine); f.setFrameShadow(QFrame.Sunken); f.setFixedWidth(4)
self.grid.addWidget(f)
w=QLineEdit('')
self.grid.addWidget(w);
@@ -145,12 +186,17 @@
def refresh(self):
val=self.getter(); axis,angle=val.toAxisAngle()
for i in (0,1,2,4):
- w=self.grid.itemAt(i).widget(); w.setText(str(axis[i] if i<3 else angle));
- if True or not w.hasFocus(): w.home(False)
+ #if True or not w.hasFocus(): w.home(False)
+ w=self.grid.itemAt(i).widget();
+ if (not w.hasFocus()):
+ w.setText(str(axis[i] if i<3 else angle));
+ w.home(False)
def update(self):
try:
x=[float((self.grid.itemAt(i).widget().text())) for i in (0,1,2,4)]
- except ValueError: self.refresh()
+ except ValueError:
+ self.refresh()
+ return
q=Quaternion(Vector3(x[0],x[1],x[2]),x[3]); q.normalize() # from axis-angle
self.trySetter(q)
def setFocus(self): self.grid.itemAt(0).widget().setFocus()
@@ -173,16 +219,24 @@
def refresh(self):
pos,ori=self.getter(); axis,angle=ori.toAxisAngle()
for i in (0,1,2,4):
- w=self.grid.itemAtPosition(1,i).widget(); w.setText(str(axis[i] if i<3 else angle));
- if True or not w.hasFocus(): w.home(False)
+ #if True or not w.hasFocus(): w.home(False)
+ w=self.grid.itemAtPosition(1,i).widget();
+ if (not w.hasFocus()):
+ w.setText(str(axis[i] if i<3 else angle));
+ w.home(False)
for i in (0,1,2):
- w=self.grid.itemAtPosition(0,i).widget(); w.setText(str(pos[i]));
- if True or not w.hasFocus(): w.home(False)
+ #if True or not w.hasFocus(): w.home(False)
+ w=self.grid.itemAtPosition(0,i).widget();
+ if (not w.hasFocus()):
+ w.setText(str(pos[i]));
+ w.home(False)
def update(self):
try:
q=[float((self.grid.itemAtPosition(1,i).widget().text())) for i in (0,1,2,4)]
v=[float((self.grid.itemAtPosition(0,i).widget().text())) for i in (0,1,2)]
- except ValueError: self.refresh()
+ except ValueError:
+ self.refresh()
+ return
qq=Quaternion(Vector3(q[0],q[1],q[2]),q[3]); qq.normalize() # from axis-angle
self.trySetter((v,qq))
def setFocus(self): self.grid.itemAtPosition(0,0).widget().setFocus()
@@ -195,6 +249,7 @@
self.rows,self.cols=rows,cols
self.idxConverter=idxConverter
self.setContentsMargins(0,0,0,0)
+ self.first=True
val=self.getter()
self.grid=QGridLayout(self); self.grid.setSpacing(0); self.grid.setMargin(0)
for row,col in itertools.product(range(self.rows),range(self.cols)):
@@ -203,12 +258,16 @@
w.textEdited.connect(self.isHot)
w.selectionChanged.connect(self.isHot)
w.editingFinished.connect(self.update)
- def refresh(self):
+ def refresh(self,force=False):
val=self.getter()
for row,col in itertools.product(range(self.rows),range(self.cols)):
w=self.grid.itemAtPosition(row,col).widget()
- w.setText(str(val[self.idxConverter(row,col)]))
- if True or not w.hasFocus: w.home(False) # make the left-most part visible, if the text is wider than the widget
+ if(self.first or force):
+ w.setText(str(val[self.idxConverter(row,col)]))
+ if (not w.hasFocus):
+ w.setText(str(val[self.idxConverter(row,col)]))
+ w.home(False) # make the left-most part visible, if the text is wider than the widget
+ self.first=False
def update(self):
try:
val=self.getter()
@@ -217,7 +276,8 @@
if w.isModified(): val[self.idxConverter(row,col)]=float(w.text())
logging.debug('setting'+str(val))
self.trySetter(val)
- except ValueError: self.refresh()
+ except ValueError:
+ self.refresh(force=True)
def setFocus(self): self.grid.itemAtPosition(0,0).widget().setFocus()
class AttrEditor_MatrixXi(AttrEditor,QFrame):
@@ -276,8 +336,8 @@
class Se3FakeType: pass
-_fundamentalEditorMap={bool:AttrEditor_Bool,str:AttrEditor_Str,int:AttrEditor_Int,float:AttrEditor_Float,Quaternion:AttrEditor_Quaternion,Vector2:AttrEditor_Vector2,Vector3:AttrEditor_Vector3,Vector6:AttrEditor_Vector6,Matrix3:AttrEditor_Matrix3,Vector6i:AttrEditor_Vector6i,Vector3i:AttrEditor_Vector3i,Vector2i:AttrEditor_Vector2i,Se3FakeType:AttrEditor_Se3}
-_fundamentalInitValues={bool:True,str:'',int:0,float:0.0,Quaternion:Quaternion((0,1,0),0.0),Vector3:Vector3.Zero,Matrix3:Matrix3.Zero,Vector6:Vector6.Zero,Vector6i:Vector6i.Zero,Vector3i:Vector3i.Zero,Vector2i:Vector2i.Zero,Vector2:Vector2.Zero,Se3FakeType:(Vector3.Zero,Quaternion((0,1,0),0.0))}
+_fundamentalEditorMap={bool:AttrEditor_Bool,str:AttrEditor_Str,int:AttrEditor_Int,float:AttrEditor_Float,complex:AttrEditor_Complex,Quaternion:AttrEditor_Quaternion,Vector2:AttrEditor_Vector2,Vector3:AttrEditor_Vector3,Vector6:AttrEditor_Vector6,Matrix3:AttrEditor_Matrix3,Vector6i:AttrEditor_Vector6i,Vector3i:AttrEditor_Vector3i,Vector2i:AttrEditor_Vector2i,Se3FakeType:AttrEditor_Se3}
+_fundamentalInitValues={bool:True,str:'',int:0,float:0.0,complex:0.0j,Quaternion:Quaternion((0,1,0),0.0),Vector3:Vector3.Zero,Matrix3:Matrix3.Zero,Vector6:Vector6.Zero,Vector6i:Vector6i.Zero,Vector3i:Vector3i.Zero,Vector2i:Vector2i.Zero,Vector2:Vector2.Zero,Se3FakeType:(Vector3.Zero,Quaternion((0,1,0),0.0))}
class SerQLabel(QLabel):
def __init__(self,parent,label,tooltip,path):
@@ -338,7 +398,7 @@
return m
vecMap={
'bool':bool,'int':int,'long':int,'Body::id_t':long,'size_t':long,
- 'Real':float,'float':float,'double':float,
+ 'Real':float,'float':float,'double':float,'complex':complex,'std::complex<Real>':complex,
'Vector6r':Vector6,'Vector6i':Vector6i,'Vector3i':Vector3i,'Vector2r':Vector2,'Vector2i':Vector2i,
'Vector3r':Vector3,'Matrix3r':Matrix3,'Se3r':Se3FakeType,
'string':str,
@@ -365,8 +425,16 @@
val=getattr(self.ser,attr) # get the value using serattr, as it might be different from what the dictionary provides (e.g. Body.blockedDOFs)
t=None
doc=getattr(self.ser.__class__,attr).__doc__;
+ # some attributes are not shown
if '|yhidden|' in doc: continue
if attr in self.ignoredAttrs: continue
+# FIXME: (Janek) Implementing Quantum Mechanics makes some DEM assumptions
+# invalid. I think that we should rethink what base class Body contains, so
+# that in QM we would not need to use this hack to hide some variables.
+# However it is great to note that only this little 'cosmetic' hack is needed
+# to make Quantum Mechanics possible in yade
+# See also: class QuantumMechanicalState, class QuantumMechanicalBody, gui/qt4/SerializableEditor.py
+ if hasattr(self.ser,"qtHide") and (attr in getattr(self.ser,"qtHide").split()): continue
if isinstance(val,list):
t=self.getListTypeFromDocstring(attr)
if not t and len(val)==0: t=(val[0].__class__,) # 1-tuple is list of the contained type
@@ -378,6 +446,7 @@
flags=int(match.group(1)) if match else 0
#logging.debug('Attr %s is of type %s'%(attr,((t[0].__name__,) if isinstance(t,tuple) else t.__name__)))
self.entries.append(self.EntryData(name=attr,T=t))
+ #print("name: "+str(attr)+"\tflag: "+str(flags))
def getDocstring(self,attr=None):
"If attr is *None*, return docstring of the Serializable itself"
doc=(getattr(self.ser.__class__,attr).__doc__ if attr else self.ser.__class__.__doc__)
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-07-02 16:18:24 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-10-07 18:41:54 +0000
@@ -51,6 +51,8 @@
vector<CellHandle> IPCells;
vector<pair<Point,Real> > imposedF;
vector<CellHandle> IFCells;
+ //Blocked cells, where pressure may be computed in undrained condition
+ vector<CellHandle> blockedCells;
//Pointers to vectors used for user defined boundary pressure
vector<Real> *pxpos, *ppval;
@@ -155,6 +157,7 @@
double averagePressure();
int getCell (double X,double Y,double Z);
double boundaryFlux(unsigned int boundaryId);
+ void setBlocked(CellHandle& cell);
vector<Real> averageFluidVelocityOnSphere(unsigned int Id_sph);
//Solver?
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-07-16 16:49:41 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-10-07 18:41:54 +0000
@@ -467,6 +467,19 @@
}
template <class Tesselation>
+void FlowBoundingSphere<Tesselation>::setBlocked(CellHandle& cell)
+{
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ cerr<<"skipping a blocked cell"<<endl;
+ if (cell->info().Pcondition=true) cell->info().p() = 0;
+ else blockedCells.push_back(cell);
+ for (int j=0; j<4; j++) {
+ (cell->info().kNorm())[j]= 0;
+ (cell->neighbor(j)->info().kNorm())[Tri.mirror_index(cell, j)]= 0;}
+}
+
+
+template <class Tesselation>
void FlowBoundingSphere<Tesselation>::computePermeability()
{
if (debugOut) cout << "----Computing_Permeability------" << endl;
@@ -486,6 +499,9 @@
for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
CellHandle& cell = *cellIt;
+ if (cell->info().blocked) {
+ setBlocked(cell);
+ cell->info().isvisited = !ref;}
Point& p1 = cell->info();
for (int j=0; j<4; j++) {
neighbourCell = cell->neighbor(j);
@@ -504,6 +520,8 @@
W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1-v0)*(v2-v0)/sqrt((v1-v0).squared_length()*(v2-v0).squared_length())),
W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
+ //FIXME: it should be possible to skip completely blocked cells, currently the problem is it segfault for undefined areas
+// if (cell->info().blocked) continue;//We don't need permeability for blocked cells, it will be set to zero anyway
pass+=1;
CVector l = p1 - p2;
@@ -846,7 +864,7 @@
int bb=-1;
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
bb++;
- if ( !cell->info().Pcondition ) {
+ if ( !cell->info().Pcondition && !cell->info().blocked) {
cell2++;
#endif
if (compressible && j==0) { previousP[bb]=cell->info().p(); }
@@ -913,7 +931,6 @@
} while ((dp_max/p_max) > tolerance /*&& j<4000*/ /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
#endif
}
-
if (debugOut) {cout << "pmax " << p_max << "; pmoy : " << p_moy << endl;
cout << "iteration " << j <<"; erreur : " << dp_max/p_max << endl;}
computedOnce=true;
=== modified file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-07-02 16:18:24 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp 2014-10-07 18:41:54 +0000
@@ -156,7 +156,7 @@
const FiniteCellsIterator cellEnd = Tri.finite_cells_end();
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
orderedCells.push_back(cell);
- if (!cell->info().Pcondition) ++ncols;}
+ if (!cell->info().Pcondition && !cell->info().blocked) ++ncols;}
// //Segfault on 14.10, and useless overall since SuiteSparse has preconditionners (including metis)
// spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
@@ -189,7 +189,7 @@
FiniteCellsIterator& cell = orderedCells[i];
///Non-ordered cells
// for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- if (!cell->info().Pcondition) {
+ if (!cell->info().Pcondition && !cell->info().blocked) {
index=cell->info().index;
if (index==0) {
T_cells[++T_index]=cell;
@@ -325,7 +325,7 @@
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
orderedCells.push_back(cell);
- if (!cell->info().Pcondition) ++ncols;
+ if (!cell->info().Pcondition && !cell->info().blocked) ++ncols;
}
//FIXME: does it really help? test by commenting this "sorting" line
spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
@@ -361,7 +361,7 @@
FiniteCellsIterator& cell = orderedCells[i];
///Non-ordered cells
// for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- if (!cell->info().Pcondition) {
+ if (!cell->info().Pcondition && !cell->info().blocked) {
if (cell->info().index==0) {
T_cells[++T_index]=cell;
cell->info().index=T_index;
@@ -403,7 +403,7 @@
FiniteCellsIterator& cell = orderedCells[i];
///Non-ordered cells
// for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- if (!cell->info().Pcondition) for (int j=0; j<4; j++) {
+ if (!cell->info().Pcondition && !cell->info().blocked) for (int j=0; j<4; j++) {
CellHandle neighbourCell = cell->neighbor(j);
if (!Tri.is_infinite(neighbourCell) && neighbourCell->info().Pcondition)
gsB[cell->info().index]+=cell->info().kNorm()[j]*neighbourCell->info().p();
=== modified file 'lib/triangulation/PeriodicFlow.hpp'
--- lib/triangulation/PeriodicFlow.hpp 2014-04-17 15:10:23 +0000
+++ lib/triangulation/PeriodicFlow.hpp 2014-10-07 18:44:58 +0000
@@ -179,6 +179,9 @@
for (VCellIterator cellIt=T[currentTes].cellHandles.begin(); cellIt!=T[currentTes].cellHandles.end(); cellIt++){
CellHandle& cell = *cellIt;
Point& p1 = cell->info();
+ if (cell->info().blocked) {
+ setBlocked(cell);}
+ if (cell->info().isGhost) {cerr<<"skipping a ghost"<<endl; continue;}
for (int j=0; j<4; j++){
neighbourCell = cell->neighbor(j);
Point& p2 = neighbourCell->info();
@@ -202,6 +205,9 @@
W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
#endif
+ //FIXME: it should be possible to skip completely blocked cells, currently the problem is it segfault for undefined areas
+ //if (cell->info().blocked) continue;//We don't need permeability for blocked cells, it will be set to zero anyway
+
pass+=1;
CVector l = p1 - p2;
distance = sqrt(l.squared_length());
@@ -219,14 +225,14 @@
const CVector& Surfk = cell->info().facetSurfaces[j];
Real area = sqrt(Surfk.squared_length());
const CVector& crossSections = cell->info().facetSphereCrossSections[j];
- Real S0=0;
- S0=checkSphereFacetOverlap(v0,v1,v2);
- if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
- if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
- //take absolute value, since in rare cases the surface can be negative (overlaping spheres)
- fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
- cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
- k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
+ Real S0=0;
+ S0=checkSphereFacetOverlap(v0,v1,v2);
+ if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
+ if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
+ //take absolute value, since in rare cases the surface can be negative (overlaping spheres)
+ fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
+ cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
+ k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
meanDistance += distance;
meanRadius += radius;
meanK += k*kFactor;
@@ -234,7 +240,7 @@
if (k<0 && debugOut) {surfneg+=1;
cout<<"__ k<0 __"<<k<<" "<<" fluidArea "<<fluidArea<<" area "<<area<<" "<<crossSections[0]<<" "<<crossSections[1]<<" "<<crossSections[2] <<" "<<W[0]->info().id()<<" "<<W[1]->info().id()<<" "<<W[2]->info().id()<<" "<<p1<<" "<<p2<<" test "<<test<<endl;}
- } else cout <<"infinite K2! surfaces will be missing (FIXME)"<<endl; k = infiniteK;
+ } else {cout <<"infinite K2! surfaces will be missing (FIXME)"<<endl; k = infiniteK;}
//Will be corrected in the next loop
(cell->info().kNorm())[j]= k*kFactor;
if (!neighbourCell->info().isGhost) (neighbourCell->info().kNorm())[Tri.mirror_index(cell, j)]= (cell->info().kNorm())[j];
=== modified file 'pkg/common/GravityEngines.cpp'
--- pkg/common/GravityEngines.cpp 2014-07-03 07:16:58 +0000
+++ pkg/common/GravityEngines.cpp 2014-09-15 10:27:11 +0000
@@ -16,7 +16,7 @@
CREATE_LOGGER(GravityEngine);
void GravityEngine::action(){
- if (warnOnce) {warnOnce=false; LOG_ERROR("GravityEngine is deprecated, consider using Newton::gravity instead (unless gravitational energy has to be tracked - not implemented with the newton attribute).")}
+ if (warnOnce) {warnOnce=false; LOG_WARN("GravityEngine is deprecated, consider using Newton::gravity instead (unless gravitational energy has to be tracked - not implemented with the newton attribute).")}
const bool trackEnergy(scene->trackEnergy);
const Real dt(scene->dt);
YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
=== modified file 'pkg/common/Grid.cpp'
--- pkg/common/Grid.cpp 2014-07-25 16:12:46 +0000
+++ pkg/common/Grid.cpp 2014-09-23 13:00:45 +0000
@@ -237,24 +237,24 @@
}
Real relPosPrev = (branch.dot(segtPrev))/(segtPrev.norm()*segtPrev.norm());
// ... and check whether the sphere projection is before the neighbours connections too.
- const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,gridNo1->ConnList[i]->getId());
if(relPosPrev<=0){ //if the sphere projection is outside both the current Connection AND this neighbouring connection, then create the interaction if the neighbour did not already do it before.
const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,gridNo1->ConnList[i]->getId());
if(intr && intr->isReal()){
shared_ptr<ScGridCoGeom> intrGeom=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
if(!(intrGeom->isDuplicate==1)){ //skip contact.
- if (isNew) {return false;}
- else {scm->isDuplicate=1;}/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/
+ if (isNew) return false;
+ else {scm->isDuplicate=1;/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/}
}
}
}
else{//the sphere projection is outside the current Connection but inside the previous neighbour. The contact has to be handled by the Prev GridConnection, not here.
if (isNew)return false;
else {
- //cout<<"The contact "<<c->id1<<"-"<<c->id2<<" HAVE to be copied and deleted NOW."<<endl ;
+// cout<<"The contact "<<c->id1<<"-"<<c->id2<<" may be copied and will be deleted now."<<endl ;
scm->isDuplicate=1;
scm->trueInt=-1;
- return true;}
+ return true;
+ }
}
}
}
@@ -279,22 +279,24 @@
shared_ptr<ScGridCoGeom> intrGeom=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
if(!(intrGeom->isDuplicate==1)){
if (isNew) return false;
- else scm->isDuplicate=1;/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/
+ else {scm->isDuplicate=1;/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/}
}
- }
+ }
}
else{//the sphere projection is outside the current Connection but inside the previous neighbour. The contact has to be handled by the Prev GridConnection, not here.
if (isNew)return false;
- else {//cout<<"The contact "<<c->id1<<"-"<<c->id2<<" HAVE to be copied and deleted NOW."<<endl ;
+ else {
+// cout<<"The contact "<<c->id1<<"-"<<c->id2<<" may be copied and will be deleted now."<<endl ;
scm->isDuplicate=1 ;
scm->trueInt=-1 ;
- return true;}
+ return true;
+ }
}
}
}
}
- else if (isNew && relPos<=0.5){
+ else if (relPos<=0.5){
if(gridNo1->ConnList.size()>1){// if the node is not an extremity of the Grid (only one connection)
for(int unsigned i=0;i<gridNo1->ConnList.size();i++){ // ... loop on all the Connections of the same Node ...
GridConnection* GC = (GridConnection*)gridNo1->ConnList[i]->shape.get();
@@ -309,11 +311,13 @@
if(relPosPrev<=0){ //the sphere projection is inside the current Connection and outide this neighbour connection.
const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,gridNo1->ConnList[i]->getId());
if( intr && intr->isReal() ){// if an ineraction exist between the sphere and the previous connection, import parameters.
- //cout<<"Copying contact geom and phys from "<<intr->id1<<"-"<<intr->id2<<" to here ("<<c->id1<<"-"<<c->id2<<")"<<endl;
scm=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
- c->geom=scm;
- c->phys=intr->phys;
- c->iterMadeReal=intr->iterMadeReal;
+ if(isNew){
+// cout<<"Copying contact geom and phys from "<<intr->id1<<"-"<<intr->id2<<" to here ("<<c->id1<<"-"<<c->id2<<")"<<endl;
+ c->geom=scm;
+ c->phys=intr->phys;
+ c->iterMadeReal=intr->iterMadeReal;
+ }
scm->trueInt=c->id2;
scm->isDuplicate=2; //command the old contact deletion.
isNew=0;
@@ -324,7 +328,7 @@
}
}
- else if (isNew && relPos>0.5){
+ else if (relPos>0.5){
if(gridNo2->ConnList.size()>1){
for(int unsigned i=0;i<gridNo2->ConnList.size();i++){
GridConnection* GC = (GridConnection*)gridNo2->ConnList[i]->shape.get();
@@ -339,11 +343,13 @@
if(relPosNext<=0){ //the sphere projection is inside the current Connection and outide this neighbour connection.
const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,gridNo2->ConnList[i]->getId());
if( intr && intr->isReal() ){// if an ineraction exist between the sphere and the previous connection, import parameters.
- //cout<<"Copying contact geom and phys from "<<intr->id1<<"-"<<intr->id2<<" to here ("<<c->id1<<"-"<<c->id2<<")"<<endl;
scm=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
- c->geom=scm;
- c->phys=intr->phys;
- c->iterMadeReal=intr->iterMadeReal;
+ if(isNew){
+// cout<<"Copying contact geom and phys from "<<intr->id1<<"-"<<intr->id2<<" to here ("<<c->id1<<"-"<<c->id2<<")"<<endl;
+ c->geom=scm;
+ c->phys=intr->phys;
+ c->iterMadeReal=intr->iterMadeReal;
+ }
scm->trueInt=c->id2;
scm->isDuplicate=2; //command the old contact deletion.
isNew=0;
@@ -410,6 +416,7 @@
if (id2!=geom->trueInt) {
//cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
if (geom->isDuplicate==2) return false;
+ return true;
}
}
Real& un=geom->penetrationDepth;
@@ -461,6 +468,7 @@
if (id2!=geom->trueInt) {
//cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
if (geom->isDuplicate==2) return false;
+ return true;
}
}
=== modified file 'pkg/common/InsertionSortCollider.cpp'
--- pkg/common/InsertionSortCollider.cpp 2014-07-03 17:20:40 +0000
+++ pkg/common/InsertionSortCollider.cpp 2014-09-15 10:27:11 +0000
@@ -58,7 +58,7 @@
void InsertionSortCollider::insertionSortParallel(VecBounds& v, InteractionContainer* interactions, Scene*, bool doCollide){
#ifdef YADE_OPENMP
assert(!periodic);
- assert(v.size==(long)v.vec.size());
+ assert(v.size==(long)v.vec.size());
if (ompThreads<=1) return insertionSort(v,interactions, scene, doCollide);
Real chunksVerlet = 4*verletDist;//is 2* the theoretical requirement?
@@ -73,16 +73,16 @@
bool changeChunks=false;
for(unsigned n=1; n<nChunks;n++) if (chunksVerlet>(v[chunks[n]].coord-v[chunks[n-1]].coord)) changeChunks=true;
if (!changeChunks) break;
- nChunks--; chunkSize = unsigned(v.size/nChunks)+1; chunks.clear();
+ nChunks--; chunkSize = unsigned(v.size/nChunks)+1; chunks.clear();
for(unsigned n=0; n<nChunks;n++) chunks.push_back(n*chunkSize); chunks.push_back(v.size);
}
static unsigned warnOnce=0;
- if (nChunks<unsigned(ompThreads) && !warnOnce++) LOG_WARN("Parallel insertion: only "<<nChunks <<" thread(s) used. The number of bodies is probably too small for allowing more threads. The contact detection should succeed but not all available threads are used.");
+ if (nChunks<unsigned(ompThreads) && !warnOnce++) LOG_WARN("Parallel insertion: only "<<nChunks <<" thread(s) used. The number of bodies is probably too small for allowing more threads, or the geometry is flat. The contact detection should succeed but not all available threads are used.");
///Define per-thread containers bufferizing the actual insertion of new interactions, since inserting is not thread-safe
std::vector<std::vector<std::pair<Body::id_t,Body::id_t> > > newInteractions;
newInteractions.resize(ompThreads,std::vector<std::pair<Body::id_t,Body::id_t> >());
- for (int kk=0; kk<ompThreads; kk++) newInteractions[kk].reserve(100);
+ for (int kk=0; kk<ompThreads; kk++) newInteractions[kk].reserve(long(chunkSize*0.3));
/// First sort, independant in each chunk
#pragma omp parallel for schedule(dynamic,1) num_threads(ompThreads>0 ? min(ompThreads,omp_get_max_threads()) : omp_get_max_threads())
@@ -103,21 +103,22 @@
v[j+1]=viInit;
}
}
-
///In the second sort, the chunks are connected consistently.
///If sorting requires to move a bound past half-chunk, the algorithm is not thread safe,
- /// if it happens we roll-back and run the 1-thread sort + send warning
+ /// if it happens we run the 1-thread sort at the end
bool parallelFailed=false;
#pragma omp parallel for schedule(dynamic,1) num_threads(ompThreads>0 ? min(ompThreads,omp_get_max_threads()) : omp_get_max_threads())
for (unsigned k=1; k<nChunks;k++) {
-
int threadNum = omp_get_thread_num();
long i=chunks[k];
- for(; v[i]<v[i-1]; i++){
+ long halfChunkStart = long(i-chunkSize*0.5);
+ long halfChunkEnd = long(i+chunkSize*0.5);
+ for(; i<halfChunkEnd; i++){
+ if (!(v[i]<v[i-1])) break; //contiguous chunks now connected consistently
const Bounds viInit=v[i]; long j=i-1; /* cache hasBB; otherwise 1% overall performance hit */ const bool viInitBB=viInit.flags.hasBB;
const bool isMin=viInit.flags.isMin;
- while(j>=chunks[k-1] && viInit<v[j]){
+ while(j>=halfChunkStart && viInit<v[j]){
v[j+1]=v[j];
if(isMin && !v[j].flags.isMin && doCollide && viInitBB && v[j].flags.hasBB && (viInit.id!=v[j].id)) {
const Body::id_t& id1 = v[j].id; const Body::id_t& id2 = viInit.id;
@@ -127,23 +128,17 @@
j--;
}
v[j+1]=viInit;
- if (j<=long(chunks[k]-chunkSize*0.5)) {
- LOG_WARN("parallel sort not guaranteed to succeed; in chunk "<<k<<" of "<<nChunks<< ", bound descending past half-chunk. Parallel colliding aborted, starting again in single thread. Consider turning collider.ompThreads=1 for not wasting CPU time.");
- parallelFailed=true;}
+ if (j<halfChunkStart) parallelFailed=true;
}
- if (i>=long(chunks[k]+chunkSize*0.5)) {
- LOG_ERROR("parallel sort not guaranteed to succeed; in chunk "<<k+1<<" of "<<nChunks<< ", bound advancing past half-chunk. Consider turning collider.ompThreads=1 for not wasting CPU time.");
- parallelFailed=true;}
+ if (i>=halfChunkEnd) parallelFailed=true;
}
- /// Check again, just to be sure...
- for (unsigned k=1; k<nChunks;k++) if (v[chunks[k]]<v[chunks[k]-1]) {
- LOG_ERROR("Parallel colliding failed, starting again in single thread. Consider turning collider.ompThreads=1 for not wasting CPU time.");
- parallelFailed=true;}
-
+ /// Now insert interactions sequentially
+ for (int n=0;n<ompThreads;n++)
+ for (size_t k=0, kend=newInteractions[n].size();k<kend;k++)
+ /*if (!interactions->found(newInteractions[n][k].first,newInteractions[n][k].second))*/ //Not needed, already checked above
+ interactions->insert(shared_ptr<Interaction>(new Interaction(newInteractions[n][k].first,newInteractions[n][k].second)));
+ /// If some bounds traversed more than a half-chunk, we complete colliding with the sequential sort
if (parallelFailed) return insertionSort(v,interactions, scene, doCollide);
-
- /// Now insert interactions sequentially
- for (int n=0;n<ompThreads;n++) for (size_t k=0, kend=newInteractions[n].size();k<kend;k++) if (!interactions->found(newInteractions[n][k].first,newInteractions[n][k].second)) interactions->insert(shared_ptr<Interaction>(new Interaction(newInteractions[n][k].first,newInteractions[n][k].second)));
#endif
}
=== modified file 'pkg/common/PeriodicEngines.hpp'
--- pkg/common/PeriodicEngines.hpp 2014-07-19 19:52:41 +0000
+++ pkg/common/PeriodicEngines.hpp 2014-09-29 15:31:49 +0000
@@ -14,6 +14,15 @@
const Real& virtNow=scene->time;
Real realNow=getClock();
const long& iterNow=scene->iter;
+
+ if((firstIterRun > 0) && (nDone==0)) {
+ if((firstIterRun > 0) && (firstIterRun == iterNow)) {
+ realLast=realNow; virtLast=virtNow; iterLast=iterNow; nDone++;
+ return true;
+ }
+ return false;
+ }
+
if (iterNow<iterLast) nDone=0;//handle O.resetTime(), all counters will be initialized again
if((nDo<0 || nDone<nDo) &&
((virtPeriod>0 && virtNow-virtLast>=virtPeriod) ||
@@ -22,6 +31,7 @@
realLast=realNow; virtLast=virtNow; iterLast=iterNow; nDone++;
return true;
}
+
if(nDone==0){
realLast=realNow; virtLast=virtNow; iterLast=iterNow; nDone++;
if(initRun) return true;
@@ -59,6 +69,7 @@
((long,iterPeriod,((void)"deactivated",0),,"Periodicity criterion using step number (deactivated if <= 0)"))
((long,nDo,((void)"deactivated",-1),,"Limit number of executions by this number (deactivated if negative)"))
((bool,initRun,false,,"Run the first time we are called as well."))
+ ((long,firstIterRun,0,,"Sets the step number, at each an engine should be executed for the first time (disabled by default)."))
((Real,virtLast,0,,"Tracks virtual time of last run |yupdate|."))
((Real,realLast,0,,"Tracks real time of last run |yupdate|."))
((long,iterLast,0,,"Tracks step number of last run |yupdate|."))
=== modified file 'pkg/dem/CapillaryPhys.hpp'
--- pkg/dem/CapillaryPhys.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/CapillaryPhys.hpp 2014-10-01 17:33:13 +0000
@@ -17,11 +17,11 @@
virtual ~CapillaryPhys() {};
- YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(CapillaryPhys,FrictPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(CapillaryPhys,FrictPhys,"Physics (of interaction) for :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`.",
((bool,meniscus,false,,"Presence of a meniscus if true"))
((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
- ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
- ((Real,vMeniscus,0.,,"Volume of the menicus"))
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc. Defined as Ugas-Uliquid, obtained from :yref:`corresponding Law2 parameter<Law2_ScGeom_CapillaryPhys_Capillarity.capillaryPressure>`"))
+ ((Real,vMeniscus,0.,,"Volume of the meniscus"))
((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
@@ -46,7 +46,7 @@
const shared_ptr<Interaction>& interaction);
FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_FrictMat_CapillaryPhys,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",,;);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_FrictMat_CapillaryPhys,IPhysFunctor, "RelationShips to use with :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`.\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",,;);
};
REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_CapillaryPhys);
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-07-18 18:18:50 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-09-17 16:40:57 +0000
@@ -87,7 +87,7 @@
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_s$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. The maximum value of moments can be defined and takes the form of rolling friction. Cohesive -type moment may also be included in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_s$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False, the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$, so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies (details can be found in [Bourrier2013]_). The maximum value of moments can be defined and takes the form of rolling friction. Cohesive -type moment may also be included in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
((bool,shear_creep,false,,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
=== modified file 'pkg/dem/HertzMindlin.hpp'
--- pkg/dem/HertzMindlin.hpp 2014-07-23 09:07:11 +0000
+++ pkg/dem/HertzMindlin.hpp 2014-10-01 17:33:13 +0000
@@ -181,8 +181,8 @@
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(MindlinCapillaryPhys,MindlinPhys,"Adds capillary physics to Mindlin's interaction physics.",
((bool,meniscus,false,,"Presence of a meniscus if true"))
((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
- ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
- ((Real,vMeniscus,0.,,"Volume of the menicus"))
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc. Defined as Ugas-Uliquid, obtained from :yref:`corresponding Law2 parameter<Law2_ScGeom_CapillaryPhys_Capillarity.capillaryPressure>`"))
+ ((Real,vMeniscus,0.,,"Volume of the meniscus"))
((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2014-10-01 17:33:13 +0000
@@ -18,7 +18,7 @@
- (french, lot of documentation) L. Scholtes, PhD thesis -> http://tel.archives-ouvertes.fr/tel-00363961/en/
- (english, less...) L. Scholtes et al. Micromechanics of granular materials with capillary effects. International Journal of Engineering Science 2009,(47)1, 64-75
-The law needs ascii files M(r=i) with i=R1/R2 to work (downloaded from https://yade-dem.org/wiki/CapillaryTriaxialTest). They contain a set of results from the resolution of the Laplace-Young equation for different configurations of the interacting geometry and must be placed in the bin directory (where yade exec file is situated) to be taken into account.
+The law needs ascii files M(r=i) with i=R1/R2 to work (downloaded from https://yade-dem.org/wiki/CapillaryTriaxialTest). They 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) and must be placed in the bin directory (where yade exec file is situated) to be taken into account.
The control parameter is the capillary pressure (or suction) Delta_u, defined as the difference between gas and liquid pressure: Delta_u = u_gas - u_liquid
Liquid bridges properties (volume V, extent over interacting grains delta1 and delta2) are computed as a result of Delta_u and the interacting geometry (spheres radii and interparticular distance)
@@ -33,10 +33,10 @@
class MeniscusParameters
{
public :
- Real V;
- Real F;
- Real delta1;
- Real delta2;
+ Real V; // adimentionnal volume of the meniscus : true volume / Rmax^3, see Annexe 1 of Scholtes2009d
+ Real F; // adimentionnal capillary force for this meniscus : true force / ( 2 * pi * Rmax * superficial tension), (30) of Annexe1 of Scholtes2009d
+ Real delta1; // angle defined Fig 2.5 Scholtes2009d
+ Real delta2; // angle defined Fig 2.5 Scholtes2009d
int index1;
int index2;
@@ -85,8 +85,8 @@
void action();
void postLoad(Law2_ScGeom_CapillaryPhys_Capillarity&);
- 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 capillary pressure (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.\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 defines as Uc=Ugas-Uliquid"))
+ 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"))
((bool,binaryFusion,true,,"If true, capillary forces are set to zero as soon as, at least, 1 overlap (menisci fusion) is detected"))
((bool,hertzOn,false,,"|yupdate| true if hertz model is used"))
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2014-10-08 05:24:38 +0000
@@ -84,25 +84,8 @@
#endif
}
-#ifdef YADE_OPENMP
-void NewtonIntegrator::ensureSync()
-{
- if (syncEnsured) return;
- YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
-// if(b->isClump()) continue;
- scene->forces.addForce(b->getId(),Vector3r(0,0,0));
- } YADE_PARALLEL_FOREACH_BODY_END();
- syncEnsured=true;
-}
-#endif
-
void NewtonIntegrator::action()
{
- #ifdef YADE_OPENMP
- //prevent https://bugs.launchpad.net/yade/+bug/923929
- ensureSync();
- #endif
-
scene->forces.sync();
bodySelected=(scene->selectedBody>=0);
if(warnNoForceReset && scene->forces.lastReset<scene->iter) LOG_WARN("O.forces last reset in step "<<scene->forces.lastReset<<", while the current step is "<<scene->iter<<". Did you forget to include ForceResetter in O.engines?");
@@ -291,7 +274,7 @@
if (ts) {
ts->densityScaling=dsc;
densityScaling=dsc;
- LOG_WARN("GlobalStiffnessTimeStepper found in O.engines and adjusted to match this setting. Revert in the the timestepper if you don't want the scaling adjusted automatically.");
+ LOG_WARN("GlobalStiffnessTimeStepper found in O.engines and adjusted to match this setting. Revert in the timestepper if you don't want the scaling adjusted automatically.");
return;
}
} LOG_WARN("GlobalStiffnessTimeStepper not found in O.engines. Density scaling will have no effect unless a scaling is specified manually for some bodies");
=== modified file 'pkg/dem/Polyhedra.cpp'
--- pkg/dem/Polyhedra.cpp 2014-07-18 18:18:50 +0000
+++ pkg/dem/Polyhedra.cpp 2014-09-22 12:28:50 +0000
@@ -13,7 +13,7 @@
#define _USE_MATH_DEFINES
-YADE_PLUGIN(/* self-contained in hpp: */ (Polyhedra) (PolyhedraGeom) (Bo1_Polyhedra_Aabb) (PolyhedraPhys) (PolyhedraMat) (Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys) (PolyhedraVolumetricLaw)
+YADE_PLUGIN(/* self-contained in hpp: */ (Polyhedra) (PolyhedraGeom) (Bo1_Polyhedra_Aabb) (PolyhedraPhys) (PolyhedraMat) (Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys) (Law2_PolyhedraGeom_PolyhedraPhys_Volumetric)
/* some code in cpp (this file): */
#ifdef YADE_OPENGL
(Gl1_Polyhedra) (Gl1_PolyhedraGeom) (Gl1_PolyhedraPhys)
@@ -220,6 +220,7 @@
}
}
+
//****************************************************************************************
/* Destructor */
@@ -395,9 +396,9 @@
//**************************************************************************************
#if 1
-Real PolyhedraVolumetricLaw::getPlasticDissipation() {return (Real) plasticDissipation;}
-void PolyhedraVolumetricLaw::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
-Real PolyhedraVolumetricLaw::elasticEnergy()
+Real Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::getPlasticDissipation() {return (Real) plasticDissipation;}
+void Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
+Real Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::elasticEnergy()
{
Real energy=0;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
@@ -413,7 +414,7 @@
//**************************************************************************************
// Apply forces on polyhedrons in collision based on geometric configuration
-bool PolyhedraVolumetricLaw::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
+bool Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
if (!I->geom) {return true;}
const shared_ptr<PolyhedraGeom>& contactGeom(YADE_PTR_DYN_CAST<PolyhedraGeom>(I->geom));
@@ -424,7 +425,8 @@
PolyhedraPhys* phys = dynamic_cast<PolyhedraPhys*>(I->phys.get());
//erase the interaction when aAbB shows separation, otherwise keep it to be able to store previous separating plane for fast detection of separation
- if (A->bound->min[0] >= B->bound->max[0] || B->bound->min[0] >= A->bound->max[0] || A->bound->min[1] >= B->bound->max[1] || B->bound->min[1] >= A->bound->max[1] || A->bound->min[2] >= B->bound->max[2] || B->bound->min[2] >= A->bound->max[2]) {
+ Vector3r shift2=scene->cell->hSize*I->cellDist.cast<Real>();
+ if (A->bound->min[0] >= B->bound->max[0]+shift2[0] || B->bound->min[0]+shift2[0] >= A->bound->max[0] || A->bound->min[1] >= B->bound->max[1]+shift2[1] || B->bound->min[1]+shift2[1] >= A->bound->max[1] || A->bound->min[2] >= B->bound->max[2]+shift2[2] || B->bound->min[2]+shift2[2] >= A->bound->max[2]) {
return false;
}
=== modified file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp 2014-07-18 18:18:50 +0000
+++ pkg/dem/Polyhedra.hpp 2014-09-22 12:28:50 +0000
@@ -240,26 +240,26 @@
//***************************************************************************
/*! Calculate physical response based on penetration configuration given by TTetraGeom. */
-class PolyhedraVolumetricLaw: public LawFunctor{
+class Law2_PolyhedraGeom_PolyhedraPhys_Volumetric: public LawFunctor{
OpenMPAccumulator<Real> plasticDissipation;
virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
Real elasticEnergy ();
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(PolyhedraVolumetricLaw,LawFunctor,"Calculate physical response of 2 :yref:`vector<Polyhedra>` in interaction, based on penetration configuration given by :yref:`PolyhedraGeom`.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_PolyhedraGeom_PolyhedraPhys_Volumetric,LawFunctor,"Calculate physical response of 2 :yref:`vector<Polyhedra>` in interaction, based on penetration configuration given by :yref:`PolyhedraGeom`.",
((Vector3r,shearForce,Vector3r::Zero(),,"Shear force from last step"))
((bool,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing"))
((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
((int,elastPotentialIx,-1,(Attr::hidden|Attr::noSave),"Index for elastic potential energy (with O.trackEnergy)"))
,,
- .def("elasticEnergy",&PolyhedraVolumetricLaw::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
- .def("plasticDissipation",&PolyhedraVolumetricLaw::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_CundallStrack::traceEnergy` is true.")
- .def("initPlasticDissipation",&PolyhedraVolumetricLaw::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
+ .def("elasticEnergy",&Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
+ .def("plasticDissipation",&Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_CundallStrack::traceEnergy` is true.")
+ .def("initPlasticDissipation",&Law2_PolyhedraGeom_PolyhedraPhys_Volumetric::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
);
FUNCTOR2D(PolyhedraGeom,PolyhedraPhys);
DECLARE_LOGGER;
};
-REGISTER_SERIALIZABLE(PolyhedraVolumetricLaw);
+REGISTER_SERIALIZABLE(Law2_PolyhedraGeom_PolyhedraPhys_Volumetric);
//***************************************************************************
=== modified file 'pkg/dem/Polyhedra_Ig2.cpp'
--- pkg/dem/Polyhedra_Ig2.cpp 2014-05-16 11:58:59 +0000
+++ pkg/dem/Polyhedra_Ig2.cpp 2014-09-22 12:28:50 +0000
@@ -34,7 +34,7 @@
//move and rotate 2nd the CGAL structure Polyhedron
rot_mat = (se32.orientation).toRotationMatrix();
- trans_vec = se32.position;
+ trans_vec = se32.position + shift2;
t_rot_trans = Transformation(rot_mat(0,0),rot_mat(0,1),rot_mat(0,2), trans_vec[0],rot_mat(1,0),rot_mat(1,1),rot_mat(1,2),trans_vec[1],rot_mat(2,0),rot_mat(2,1),rot_mat(2,2),trans_vec[2],1.);
Polyhedron PB = B->GetPolyhedron();
std::transform( PB.points_begin(), PB.points_end(), PB.points_begin(), t_rot_trans);
@@ -57,18 +57,23 @@
//find intersection Polyhedra
Polyhedron Int;
- Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position), bang->sep_plane);
+ Int = Polyhedron_Polyhedron_intersection(PA,PB,ToCGALPoint(bang->contactPoint),ToCGALPoint(se31.position),ToCGALPoint(se32.position+shift2), bang->sep_plane);
//volume and centroid of intersection
double volume;
Vector3r centroid;
P_volume_centroid(Int, &volume, ¢roid);
- if(isnan(volume) || volume<=1E-25 || volume > min(A->GetVolume(),B->GetVolume())) {bang->equivalentPenetrationDepth=0; return true;}
+ if(isnan(volume) || volume<=1E-25 || volume > min(A->GetVolume(),B->GetVolume())) {
+ bang->equivalentPenetrationDepth=0;
+ bang->penetrationVolume=min(A->GetVolume(),B->GetVolume());
+ bang->normal = (A->GetVolume()>B->GetVolume() ? 1 : -1)*(se32.position+shift2-se31.position);
+ return true;
+ }
if ( (!Is_inside_Polyhedron(PA, ToCGALPoint(centroid))) or (!Is_inside_Polyhedron(PB, ToCGALPoint(centroid)))) {bang->equivalentPenetrationDepth=0; return true;}
//find normal direction
Vector3r normal = FindNormal(Int, PA, PB);
- if((se32.position-centroid).dot(normal)<0) normal*=-1;
+ if((se32.position+shift2-centroid).dot(normal)<0) normal*=-1;
//calculate area of projection of Intersection into the normal plane
//double area = CalculateProjectionArea(Int, ToCGALVector(normal));
@@ -93,7 +98,6 @@
PrintPolyhedron2File(Int,fin);
fclose(fin);
*/
-
return true;
}
=== modified file 'pkg/dem/Polyhedra_Ig2.hpp'
--- pkg/dem/Polyhedra_Ig2.hpp 2013-10-16 15:24:49 +0000
+++ pkg/dem/Polyhedra_Ig2.hpp 2014-09-22 12:28:50 +0000
@@ -13,6 +13,7 @@
public:
virtual ~Ig2_Polyhedra_Polyhedra_PolyhedraGeom(){};
virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c) { return go(cm1,cm2,state2,state1,-shift2,force,c); }
FUNCTOR2D(Polyhedra,Polyhedra);
DEFINE_FUNCTOR_ORDER_2D(Polyhedra,Polyhedra);
YADE_CLASS_BASE_DOC(Ig2_Polyhedra_Polyhedra_PolyhedraGeom,IGeomFunctor,"Create/update geometry of collision between 2 Polyhedras");
=== modified file 'pkg/dem/SpherePack.cpp'
--- pkg/dem/SpherePack.cpp 2014-09-03 15:59:54 +0000
+++ pkg/dem/SpherePack.cpp 2014-09-10 20:45:19 +0000
@@ -276,17 +276,20 @@
long SpherePack::particleSD2(const vector<Real>& radii, const vector<Real>& passing, int numSph, bool periodic, Real cloudPorosity, int seed){
//deprecated (https://bugs.launchpad.net/yade/+bug/1024443)
LOG_ERROR("particleSD2() has been removed. Please use makeCloud() instead.");
+ return 1;
};
// Discrete particle size distribution
long SpherePack::particleSD(Vector3r mn, Vector3r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount, int seed){
//deprecated (https://bugs.launchpad.net/yade/+bug/1024443)
LOG_ERROR("particleSD() has been removed. Please use makeCloud() instead.");
+ return 1;
}
long SpherePack::particleSD_2d(Vector2r mn, Vector2r mx, Real rMean, bool periodic, string name, int numSph, const vector<Real>& radii, const vector<Real>& passing, bool passingIsNotPercentageButCount, int seed){
//deprecated (https://bugs.launchpad.net/yade/+bug/1024443)
LOG_ERROR("particleSD_2d() has been removed. Please use makeCloud() instead.");
+ return 1;
}
long SpherePack::makeClumpCloud(const Vector3r& mn, const Vector3r& mx, const vector<shared_ptr<SpherePack> >& _clumps, bool periodic, int num, int seed){
=== modified file 'pkg/dem/Tetra.cpp'
--- pkg/dem/Tetra.cpp 2014-07-18 18:18:50 +0000
+++ pkg/dem/Tetra.cpp 2014-09-22 12:28:50 +0000
@@ -513,7 +513,7 @@
for (int i=0; i<4; i++) {
vTemp = se31.position + se31.orientation*shape1->v[i];
p1[i] = Point(vTemp[0],vTemp[1],vTemp[2]);
- vTemp = se32.position + se32.orientation*shape2->v[i];
+ vTemp = se32.position + se32.orientation*shape2->v[i] + shift2;
p2[i] = Point(vTemp[0],vTemp[1],vTemp[2]);
}
=== modified file 'pkg/dem/TriaxialStressController.cpp'
--- pkg/dem/TriaxialStressController.cpp 2014-08-13 13:31:41 +0000
+++ pkg/dem/TriaxialStressController.cpp 2014-10-08 12:04:24 +0000
@@ -17,6 +17,11 @@
#include<yade/pkg/dem/Shop.hpp>
#include<yade/core/Clump.hpp>
+#ifdef FLOW_ENGINE
+//#include<yade/pkg/pfv/FlowEngine.hpp>
+#include "FlowEngine_FlowEngineT.hpp"
+#endif
+
CREATE_LOGGER(TriaxialStressController);
YADE_PLUGIN((TriaxialStressController));
@@ -32,8 +37,17 @@
);
}
-void TriaxialStressController::updateStiffness ()
-{
+void TriaxialStressController::updateStiffness() {
+ Real fluidStiffness = 0.;
+ #ifdef FLOW_ENGINE
+ FOREACH(const shared_ptr<Engine> e, Omega::instance().getScene()->engines) {
+ if (e->getClassName() == "FlowEngine") {
+ TemplateFlowEngine_FlowEngineT<FlowCellInfo_FlowEngineT,FlowVertexInfo_FlowEngineT>* flow =
+ dynamic_cast<TemplateFlowEngine_FlowEngineT<FlowCellInfo_FlowEngineT,FlowVertexInfo_FlowEngineT>*>(e.get());
+ if ( (flow->fluidBulkModulus > 0) && (!(flow->dead)) ) fluidStiffness = flow->fluidBulkModulus/porosity;
+ }
+ }
+ #endif
for (int i=0; i<6; ++i) stiffness[i] = 0;
InteractionContainer::iterator ii = scene->interactions->begin();
InteractionContainer::iterator iiEnd = scene->interactions->end();
@@ -46,12 +60,19 @@
int id1 = contact->getId1(), id2 = contact->getId2();
for (int index=0; index<6; ++index) if ( wall_id[index]==id1 || wall_id[index]==id2 )
{
- FrictPhys* currentContactPhysics =
- static_cast<FrictPhys*> ( contact->phys.get() );
- stiffness[index] += currentContactPhysics->kn;
+ FrictPhys* currentContactPhysics = static_cast<FrictPhys*> ( contact->phys.get() );
+ stiffness[index] += currentContactPhysics->kn;
}
}
}
+ if (fluidStiffness > 0) {
+ stiffness[0] += fluidStiffness*width*depth/height;
+ stiffness[1] += fluidStiffness*width*depth/height;
+ stiffness[2] += fluidStiffness*height*depth/width;
+ stiffness[3] += fluidStiffness*height*depth/width;
+ stiffness[4] += fluidStiffness*width*height/depth;
+ stiffness[5] += fluidStiffness*width*height/depth;
+ }
}
void TriaxialStressController::controlExternalStress(int wall, Vector3r resultantForce, State* p, Real wall_max_vel)
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2014-08-20 13:03:23 +0000
+++ pkg/dem/VTKRecorder.cpp 2014-10-07 13:13:28 +0000
@@ -662,20 +662,6 @@
const Box* box = dynamic_cast<Box*>(b->shape.get());
if (box){
- if(recActive[REC_FORCE]){
- scene->forces.sync();
- const Vector3r& f = scene->forces.getForce(b->getId());
- const Vector3r& t = scene->forces.getTorque(b->getId());
- Real ff[3] = { (Real) f[0], (Real) f[1], (Real) f[2] };
- Real tt[3] = { (Real) t[0], (Real) t[1], (Real) t[2] };
- Real fn = f.norm();
- Real tn = t.norm();
- boxesForceVec->InsertNextTupleValue(ff);
- boxesForceLen->InsertNextValue(fn);
- boxesTorqueVec->InsertNextTupleValue(tt);
- boxesTorqueLen->InsertNextValue(tn);
- }
-
Vector3r pos(scene->isPeriodic ? scene->cell->wrapShearedPt(b->state->pos) : b->state->pos);
Vector3r ext(box->extents);
vtkSmartPointer<vtkQuad> boxes = vtkSmartPointer<vtkQuad>::New();
=== modified file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 2014-07-21 09:09:14 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-09-11 17:15:11 +0000
@@ -142,12 +142,15 @@
Vector3r torque1 = Vector3r::Zero();
Vector3r torque2 = Vector3r::Zero();
- computeForceTorqueViscEl(_geom, _phys, I, force, torque1, torque2);
-
- addForce (id1,-force,scene);
- addForce (id2, force,scene);
- addTorque(id1, torque1,scene);
- addTorque(id2, torque2,scene);
+ if (computeForceTorqueViscEl(_geom, _phys, I, force, torque1, torque2)) {
+ addForce (id1,-force,scene);
+ addForce (id2, force,scene);
+ addTorque(id1, torque1,scene);
+ addTorque(id2, torque2,scene);
+ return true;
+ } else {
+ return false;
+ }
}
return true;
}
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-06-20 15:28:21 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-10-07 18:41:54 +0000
@@ -1,4 +1,5 @@
#pragma once
+
/// Frequently used:
typedef CGT::CVector CVector;
typedef CGT::Point Point;
@@ -9,6 +10,18 @@
inline Vector3r makeVector3r ( const Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
inline Vector3r makeVector3r ( const CVector& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
+/// The following macros can be used to expose CellInfo members.
+/// The syntax is CELL_SCALAR_GETTER(double,.p(),pressure), note the "." before member name, data members would be without the "()"
+#define CELL_SCALAR_GETTER(type, param, getterName) \
+type getterName(unsigned int id){\
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return 0;}\
+ return solver->T[solver->currentTes].cellHandles[id]->info()param;}
+
+#define CELL_VECTOR_GETTER(param, getterName) \
+Vector3r getterName(unsigned int id){\
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return Vector3r(0,0,0);}\
+ return makeVector3r(solver->T[solver->currentTes].cellHandles[id]->info()param);}
+
#ifdef LINSOLV
#define DEFAULTSOLVER CGT::FlowBoundingSphereLinSolv<_Tesselation>
#else
=== modified file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in 2014-09-03 15:59:54 +0000
+++ pkg/pfv/FlowEngine.hpp.in 2014-10-07 18:41:54 +0000
@@ -175,12 +175,21 @@
double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}
unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}
+ CELL_VECTOR_GETTER(,cellCenter)
+ CELL_VECTOR_GETTER(.averageCellVelocity,cellVelocity)
+ CELL_SCALAR_GETTER(double,.p(),pressure)
boost::python::list getVertices(unsigned int id){
boost::python::list ids;
- if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->vertex(i)->info().id());
return ids;
}
+ void blockCell(unsigned int id, bool blockPressure){
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return;}
+ solver->T[solver->currentTes].cellHandles[id]->info().blocked=true;
+ solver->T[solver->currentTes].cellHandles[id]->info().Pcondition=blockPressure;
+ }
+
double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
double averagePressure(){return solver->averagePressure();}
@@ -213,7 +222,7 @@
if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes();/* LOG_WARN("Computing all volumes now, as you did not request it explicitely.");*/}
return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
- YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine_@TEMPLATE_FLOW_NAME@,@TEMPLATE_FLOW_NAME@,PartialEngine,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. See the example script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n.. note::Multi-threading seems to work fine for Cholesky decomposition, but it fails for the solve phase in which -j1 is the fastest, here we specify thread numbers independently using :yref:`FlowEngine::numFactorizeThreads` and :yref:`FlowEngine::numSolveThreads`. These multhreading settings are only impacting the behaviour of openblas library and are relatively independant of :yref:`FlowEngine::multithread`. However, the settings have to be globally consistent. For instance, :yref:`multithread<FlowEngine::multithread>` =True with yref:`numFactorizeThreads<FlowEngine::numFactorizeThreads>` = yref:`numSolveThreads<FlowEngine::numSolveThreads>` = 4 implies that openblas will mobilize 8 processors at some point. If the system does not have so many procs. it will hurt performance.",
+ YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine_@TEMPLATE_FLOW_NAME@,@TEMPLATE_FLOW_NAME@,PartialEngine,"A generic engine from wich more specialized engines can inherit. It is defined for the sole purpose of inserting the right data classes CellInfo and VertexInfo in the triangulation, and it should not be used directly. Instead, look for specialized engines, e.g. :yref:`FlowEngine`, :yref:`PeriodicFlowEngine`, or :yref:`DFNFlowEngine`.",
((bool,isActivated,true,,"Activates Flow Engine"))
((bool,first,true,,"Controls the initialization/update phases"))
((bool,doInterpolate,false,,"Force the interpolation of cell's info while remeshing. By default, interpolation would be done only for compressible fluids. It can be forced with this flag."))
@@ -274,6 +283,7 @@
#endif
((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryXPos`"))
((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryPressure`"))
+ ((string,blockHook,"",,"Python command to be run after triangulation to define blocked cells (see also :yref:`TemplateFlowEngine_@TEMPLATE_FLOW_NAME@.blockCell`)"))
,
/*deprec*/
((meanK_opt,clampKValues,"the name changed"))
@@ -301,7 +311,6 @@
.def("getConstrictionsFull",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getConstrictionsFull,(boost::python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
.def("edgeSize",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::edgeSize,"Return the number of interactions.")
.def("OSI",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::OSI,"Return the number of interactions only between spheres.")
-
.def("saveVtk",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
.def("avFlVelOnSph",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")
.def("fluidForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::fluidForce,(boost::python::arg("idSph")),"Return the fluid force on sphere idSph.")
@@ -322,7 +331,9 @@
.def("updateBCs",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
.def("emulateAction",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
.def("getCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCell,(boost::python::arg("pos")),"get id of the cell containing (X,Y,Z).")
+ .def("getCellCenter",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cellCenter,(boost::python::arg("id")),"get voronoi center of cell 'id'.")
.def("nCells",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::nCells,"get the total number of finite cells in the triangulation.")
+ .def("blockCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::blockCell,(boost::python::arg("id"),boost::python::arg("blockPressure")),"block cell 'id'. The cell will be excluded from the fluid flow problem and the conductivity of all incident facets will be null. If blockPressure=False, deformation is reflected in the pressure, else it is constantly 0.")
.def("getVertices",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVertices,(boost::python::arg("id")),"get the vertices of a cell")
#ifdef LINSOLV
.def("exportMatrix",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::exportMatrix,(boost::python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
@@ -336,7 +347,7 @@
.def("averageVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageVelocity,"measure the mean velocity in the period")
)
};
-// Definition of functions in a separate file for clarity
+
class FlowCellInfo_@TEMPLATE_FLOW_NAME@ : public CGT::SimpleCellInfo {
public:
@@ -344,6 +355,7 @@
unsigned int index;
int volumeSign;
bool Pcondition;
+ bool blocked;//exclude cell from the fluid domain
Real invVoidV;
Real t;
int fict;
@@ -377,8 +389,9 @@
rayHydr.resize(4, 0);
invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;
isFictious=false; Pcondition = false; isGhost = false;
- isvisited = false;
+ isvisited = false;
isGhost=false;
+ blocked=false;
}
bool isGhost;
double invSumK;
@@ -414,4 +427,5 @@
inline const CVector ghostShift (void) const {return CGAL::NULL_VECTOR;}
};
+// Implementation of functions in this separate file for clarity
#include "FlowEngine_@TEMPLATE_FLOW_NAME@.ipp"
=== removed file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 2014-06-20 15:28:21 +0000
+++ pkg/pfv/FlowEngine.ipp 1970-01-01 00:00:00 +0000
@@ -1,1 +0,0 @@
-
=== modified file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in 2014-08-01 11:16:21 +0000
+++ pkg/pfv/FlowEngine.ipp.in 2014-10-07 18:41:54 +0000
@@ -12,6 +12,8 @@
#ifdef LINSOLV
#include <cholmod.h>
#endif
+//For pyRunString
+#include<yade/lib/pyutil/gil.hpp>
template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine_@TEMPLATE_FLOW_NAME@() {}
@@ -211,6 +213,7 @@
flow.fluidBulkModulus = fluidBulkModulus;
// flow.tesselation().Clear();
flow.tesselation().maxId=-1;
+ flow.blockedCells.clear();
flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
}
@@ -255,6 +258,7 @@
flow.tesselation().cellHandles.push_back(cell);
cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
flow.displayStatistics ();
+ if(!blockHook.empty()){ LOG_INFO("Running blockHook: "<<blockHook); pyRunString(blockHook); }
flow.computePermeability();
//This virtual function does nothing yet, derived class may overload it to make permeability different (see DFN engine)
trickPermeability();
=== modified file 'py/post2d.py'
--- py/post2d.py 2013-03-08 21:26:47 +0000
+++ py/post2d.py 2014-09-15 15:11:41 +0000
@@ -180,7 +180,12 @@
Scalar fields contain 'val' (value from *extractor*), vector fields have 'valX' and 'valY' (2 components returned by the *extractor*).
"""
- from miniEigen import Vector3
+
+ try:
+ from miniEigen import Vector3
+ except ImportError:
+ from minieigen import Vector3
+
xx,yy,dd1,dd2,rr=[],[],[],[],[]
nDim=0
objects=O.interactions if intr else O.bodies
=== modified file 'py/utils.py'
--- py/utils.py 2014-08-21 14:31:08 +0000
+++ py/utils.py 2014-09-17 00:20:31 +0000
@@ -936,7 +936,9 @@
if ((2*b.shape.radius) > maxD) : maxD = 2*b.shape.radius
if (((2*b.shape.radius) < minD) or (minD==0.0)): minD = 2*b.shape.radius
- if (minD==maxD): return false #All particles are having the same size
+ if (minD==maxD):
+ print 'Monodisperse packing with diameter =', minD,'. Not computing psd'
+ return False #All particles are having the same size
binsSizes = numpy.linspace(minD, maxD, bins+1)
=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp 2014-07-14 20:08:34 +0000
+++ py/wrapper/yadeWrapper.cpp 2014-09-26 21:38:22 +0000
@@ -830,7 +830,7 @@
.add_property("dynDtAvailable",&pyOmega::dynDtAvailable_get,"Whether a :yref:`TimeStepper` is amongst :yref:`O.engines<Omega.engines>`, activated or not.")
.def("load",&pyOmega::load,(py::arg("file"),py::arg("quiet")=false),"Load simulation from file. The file should be :yref:`saved<Omega.save>` in the same version of Yade, otherwise compatibility is not guaranteed.")
.def("reload",&pyOmega::reload,(py::arg("quiet")=false),"Reload current simulation")
- .def("save",&pyOmega::save,(py::arg("file"),py::arg("quiet")=false),"Save current simulation to file (should be .xml or .xml.bz2). The file should be :yref:`loaded<Omega.load>` in the same version of Yade, otherwise compatibility is not guaranteed.")
+ .def("save",&pyOmega::save,(py::arg("file"),py::arg("quiet")=false),"Save current simulation to file (should be .xml or .xml.bz2 or .yade or .yade.gz). .xml files are bigger than .yade, but can be more or less easily (due to their size) opened and edited, e.g. with text editors. .bz2 and .gz correspond both to compressed versions. All saved files should be :yref:`loaded<Omega.load>` in the same version of Yade, otherwise compatibility is not guaranteed.")
.def("loadTmp",&pyOmega::loadTmp,(py::arg("mark")="",py::arg("quiet")=false),"Load simulation previously stored in memory by saveTmp. *mark* optionally distinguishes multiple saved simulations")
.def("saveTmp",&pyOmega::saveTmp,(py::arg("mark")="",py::arg("quiet")=false),"Save simulation to memory (disappears at shutdown), can be loaded later with loadTmp. *mark* optionally distinguishes different memory-saved simulations.")
.def("lsTmp",&pyOmega::lsTmp,"Return list of all memory-saved simulations.")
