yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11386
[Branch ~yade-pkg/yade/git-trunk] Rev 3382: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Bruno Chareyre (bruno-chareyre)
Christian Jakob (jakob-ifgt)
Jan Stránský (honzik)...
------------------------------------------------------------
revno: 3382 [merge]
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Tue 2013-07-30 19:03:30 +0200
message:
Merge github.com:yade/trunk into chaoUnsat
Conflicts:
lib/triangulation/Network.hpp
lib/triangulation/Network.ipp
pkg/dem/FlowEngine.hpp
removed:
pkg/dem/CohesiveStateRPMRecorder.cpp
pkg/dem/CohesiveStateRPMRecorder.hpp
pkg/dem/ParticleSizeDistrbutionRPMRecorder.cpp
pkg/dem/ParticleSizeDistrbutionRPMRecorder.hpp
pkg/dem/RockPM.cpp
pkg/dem/RockPM.hpp
added:
examples/FluidCouplingPFV/
examples/FluidCouplingPFV/oedometer.py
examples/LudingPM.py
examples/PIDController.py
pkg/dem/LudingPM.cpp
pkg/dem/LudingPM.hpp
modified:
CMakeLists.txt
core/Clump.cpp
doc/references.bib
doc/yade-articles.bib
doc/yade-conferences.bib
examples/clumps/addToClump-example.py
examples/clumps/releaseFromClump-example.py
examples/test/vtk-exporter/vtkExporter.py
examples/triax-tutorial/script-session1.py
lib/base/Math.hpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/Network.hpp
lib/triangulation/Network.ipp
lib/triangulation/def_types.h
pkg/common/KinematicEngines.cpp
pkg/common/KinematicEngines.hpp
pkg/dem/BubbleMat.hpp
pkg/dem/FlowEngine.cpp
pkg/dem/FlowEngine.hpp
pkg/dem/PeriIsoCompressor.cpp
pkg/dem/PeriIsoCompressor.hpp
pkg/dem/TriaxialStressController.cpp
pkg/dem/TriaxialStressController.hpp
pkg/dem/VTKRecorder.cpp
py/export.py
py/mathWrap/miniEigen.cpp
py/tests/clump.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 'CMakeLists.txt'
--- CMakeLists.txt 2013-07-02 15:15:33 +0000
+++ CMakeLists.txt 2013-07-18 12:40:45 +0000
@@ -329,18 +329,13 @@
IF (CHUNKSIZE)
INCLUDE(CombineSources)
-
COMBINE_SOURCES(${CMAKE_BINARY_DIR}/core "${SRC_CORE}" ${CHUNKSIZE})
FILE(GLOB SRC_CORE_COMBINED "${CMAKE_BINARY_DIR}/core.*.cpp")
- ADD_LIBRARY(core SHARED ${SRC_CORE_COMBINED})
-
COMBINE_SOURCES(${CMAKE_BINARY_DIR}/pkg "${SRC_PKG}" ${CHUNKSIZE})
FILE(GLOB SRC_PKG_COMBINED "${CMAKE_BINARY_DIR}/pkg.*.cpp")
- ADD_LIBRARY(plugins SHARED ${SRC_PKG_COMBINED})
-
COMBINE_SOURCES(${CMAKE_BINARY_DIR}/lib "${SRC_LIB}" ${CHUNKSIZE})
FILE(GLOB SRC_LIB_COMBINED "${CMAKE_BINARY_DIR}/lib.*.cpp")
- ADD_LIBRARY(support SHARED ${SRC_LIB_COMBINED})
+ ADD_LIBRARY(yade SHARED ${SRC_LIB_COMBINED} ${SRC_CORE_COMBINED} ${SRC_PKG_COMBINED})
ELSE (CHUNKSIZE)
ADD_LIBRARY(yade SHARED ${SRC_CORE} ${SRC_PKG} ${SRC_LIB})
ENDIF (CHUNKSIZE)
=== modified file 'core/Clump.cpp'
--- core/Clump.cpp 2013-07-04 14:20:19 +0000
+++ core/Clump.cpp 2013-07-05 07:33:47 +0000
@@ -204,8 +204,8 @@
LOG_TRACE("R_g2c=\n"<<R_g2c);
// set quaternion from rotation matrix
state->ori=Quaternionr(R_g2c); state->ori.normalize();
- state->inertia=decomposed.eigenvalues();
- state->mass=M*dens;
+ state->inertia=dens*decomposed.eigenvalues();
+ state->mass=dens*M;
// TODO: these might be calculated from members... but complicated... - someone needs that?!
state->vel=state->angVel=Vector3r::Zero();
=== modified file 'doc/references.bib'
--- doc/references.bib 2013-06-25 04:00:41 +0000
+++ doc/references.bib 2013-07-08 08:52:03 +0000
@@ -599,3 +599,19 @@
volume = "49",
year = "2013"
}
+
+@article{Luding2008,
+ year={2008},
+ issn={1434-5021},
+ journal={Granular Matter},
+ volume={10},
+ number={4},
+ doi={10.1007/s10035-008-0099-x},
+ title={Cohesive, frictional powders: contact models for tension},
+ url={http://dx.doi.org/10.1007/s10035-008-0099-x},
+ publisher={Springer-Verlag},
+ keywords={Granular materials; Molecular dynamics (MD) and discrete elementmodel (DEM) force-laws; Friction; Rolling- and torsion-resistance; Adhesion; Plastic deformation},
+ author={Luding, Stefan},
+ pages={235-246},
+ language={English}
+}
=== modified file 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2013-06-25 04:00:41 +0000
+++ doc/yade-articles.bib 2013-07-29 17:37:33 +0000
@@ -1,13 +1,310 @@
-@article{Smilauer2006,
- title = {The splendors and miseries of Yade design},
- author = {Václav Šmilauer},
+
+@article{Boon2012a,
+ title = {A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method},
+ journal = {Computers and Geotechnics},
+ volume = {44},
+ number = {0},
+ pages = {73 - 82},
+ year = {2012},
+ doi = {10.1016/j.compgeo.2012.03.012},
+ url = {http://www.sciencedirect.com/science/article/pii/S0266352X12000535},
+ author = {Boon,C.W. and Houlsby, G.T. and Utili, S.}
+}
+
+@article{Boon2012b,
+ title = {A new contact detection algorithm for three-dimensional non-spherical particles},
+ journal = {Powder Technology},
+ year = {2013},
+ doi = {10.1016/j.powtec.2012.12.040},
+ url = {http://www.sciencedirect.com/science/article/pii/S003259101200839X},
+ author = {Boon,C.W. and Houlsby, G.T. and Utili, S.}
+}
+
+@Article{Catalano2013a,
+ title = {Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects},
+ author = {Catalano, E. and Chareyre, B. and Barthélémy, E.},
+ journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
+ year = {2013},
+ doi = {10.1002/nag.2198},
+ note = {in press, arxiv version available},
+ url = {http://arxiv.org/pdf/1304.4895.pdf}
+}
+
+@article {Chareyre2012a,
+ author = {Chareyre, B. and Cortis, A. and Catalano, E. and Barthélemy, E.},
+ title = {Pore-Scale Modeling of Viscous Flow and Induced Forces in Dense Sphere Packings},
+ journal = {Transport in Porous Media},
+ publisher = {Springer Netherlands},
+ volume={92},
+ number={2},
+ year={2012},
+ pages = {473-493},
+ url = {http://dx.doi.org/10.1007/s11242-011-9915-6},
+ doi = {10.1007/s11242-011-9915-6}
+}
+
+@article{Chen2007,
+ title = {Prediction/Verification of Particle Motion in One Dimension with the Discrete-Element Method},
+ author = {Chen, F. and Drumm, E. C. and Guiochon, G.},
+ doi = {10.1061/(ASCE)1532-3641(2007)7:5(344)},
+ journal = {International Journal of Geomechanics, ASCE},
+ pages = {344--352},
+ volume = {7},
+ number = {5},
+ year = {2007}
+}
+
+@article{Chen2011a,
+ title={Coupled discrete element and finite volume solution of two classical soil mechanics problems},
+ author={Chen, F. and Drumm, E. and Guiochon G.},
+ journal={Computers and Geotechnics},
+ year={2011},
+ doi = {10.1016/j.compgeo.2011.03.009},
+ url={http://www.sciencedirect.com/science/article/pii/S0266352X11000504},
+ publisher={Elsevier}
+}
+
+@article{Chen2012,
+ author = {Chen, Jingsong and Huang, Baoshan and Chen, Feng and Shu, Xiang},
+ title = {Application of discrete element method to Superpave gyratory compaction},
+ journal = {Road Materials and Pavement Design},
+ volume = {13},
+ number = {3},
+ pages = {480-500},
+ year = {2012},
+ doi = {10.1080/14680629.2012.694160},
+ URL = {http://www.tandfonline.com/doi/abs/10.1080/14680629.2012.694160}
+}
+
+@article{Dang2010a,
+ author = {Dang, H. K. and Meguid, M. A.},
+ title = {Algorithm to Generate a Discrete Element Specimen with Predefined Properties},
+ publisher = {ASCE},
+ year = {2010},
+ journal = {International Journal of Geomechanics},
+ volume = {10},
+ number = {2},
+ pages = {85-91},
+ keywords = {Discrete elements; Geotechnical engineering; Grain size; Porosity; Algorithms; Granular materials},
+ doi = {10.1061/(ASCE)GM.1943-5622.0000028}
+}
+
+@article{Dang2010b,
+ title = {Evaluating the performance of an explicit dynamic relaxation technique in analyzing non-linear geotechnical engineering problems},
+ journal = {Computers and Geotechnics},
+ volume = {37},
+ number = {1-2},
+ pages = {125 - 131},
+ year = {2010},
+ doi = {10.1016/j.compgeo.2009.08.004},
+ author = {Dang, H. K. and Meguid, M. A.}
+}
+
+@article{Donze2008,
+ title={Impacts on cohesive frictional geomaterials},
+ author={Donzé, F.V.},
+ journal={European Journal of Environmental and Civil Engineering},
+ volume={12},
+ number={7-8},
+ pages={967--985},
+ year={2008}
+}
+
+@article{Duriez2011a,
+ author = {Duriez,J. and Darve, F. and Donzé, F.V.},
+ title = {A discrete modeling-based constitutive relation for infilled rock joints},
+ year = {2011},
+ journal = {International Journal of Rock Mechanics & Mining Sciences},
+ volume = {48},
+ number = {3},
+ pages = {458--468},
+ doi = {10.1016/j.ijrmms.2010.09.008}
+}
+
+@article{Duriez2011b,
+ title={Incrementally non-linear plasticity applied to rock joint modelling},
+ author={Duriez, J. and Darve, F. and Donzé, F.V.},
+ journal={International Journal for Numerical and Analytical Methods in Geomechanics},
+ year={2011},
+ doi = {10.1002/nag.1105}
+}
+
+@Article{Duriez2011c,
+ author = {Duriez,J. and Darve, F. and Donzé, F.V.},
+ title = {A discrete modeling-based constitutive relation for infilled rock joints},
+ year = {2011},
+ journal = {International Journal of Rock Mechanics \& Mining Sciences},
+ volume = {48},
+ number = {3},
+ pages = {458--468},
+ doi = {10.1016/j.ijrmms.2010.09.008}
+}
+
+@article{Duriez2013,
+ author = {Duriez, J. and Darve, F. and Donzé, F.V.},
+ title = {Incrementally non-linear plasticity applied to rock joint modelling},
+ journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
+ volume = {37},
+ number = {5},
+ url = {http://dx.doi.org/10.1002/nag.1105},
+ doi = {10.1002/nag.1105},
+ pages = {453--477},
+ year = {2013}
+}
+
+@article{Epifantsev2012,
+ title={Proizvodstvo kuskovogo torfa, ekstrudirovanie, forma zakhodnoi i kalibriruyushchei chasti fil'ery matritsy, metod diskretnykh elementov [RUS]},
+ author={Epifantsev, K. and Mikhailov, A. and Gladky, A.},
+ journal={Mining informational and analytical bulletin (scientific and technical journal)},
+ year={2012},
+ publisher={The joint-stock company Gornaya Kniga},
+ number = {3},
+ pages = {212-219},
+ issn = {0236-1493},
+}
+
+@article{Favier2009a,
+ title={Predicting the drag coefficient of a granular flow using the discrete element method},
+ author={Favier, L. and Daudon, D. and Donzé, F.V. and Mazars, J.},
+ journal={Journal of Statistical Mechanics: Theory and Experiment},
+ volume={2009},
+ pages={P06012},
+ year={2009}
+}
+
+@article{Favier2012,
+ title={Rigid obstacle impacted by a supercritical cohesive granular flow using a 3D discrete element model},
+ author={Favier, L. and Daudon, D. and Donzé, F.V.},
+ journal={Cold Regions Science and Technology},
+ year={2012},
+ volume = {85},
+ pages = {232--241},
+ url = {http://dx.doi.org/10.1016/j.coldregions.2012.09.010}
+}
+
+@article{Grujicic2013,
+ title={Discrete Element Modeling and Analysis of Structural Collapse/Survivability of a Building Subjected to Improvised Explosive Device (IED) Attack},
+ author={Grujicic, M and Snipes, JS and Ramaswami, S and Yavari, R},
+ journal={Advances in Materials Science and Applications},
+ year={2013},
+ volume={2},
+ number={1},
+ pages={9--24}
+}
+
+@article{Gusenbauer2012,
+ title={Self-organizing magnetic beads for biomedical applications},
+ author={Gusenbauer, M. and Kovacs, A. and Reichel, F. and Exl, L. and Bance, S. and {\"O}zelt, H. and Schrefl, T.},
+ journal={Journal of Magnetism and Magnetic Materials},
+ volume={324},
+ number={6},
+ pages={977--982},
+ year={2012},
+ publisher={Elsevier}
+}
+
+@article{Harthong2009,
+ author = {Harthong, B. and Jerier, J.F. and Doremus, P. and Imbault, D. and Donzé, F.V.},
+ doi = {10.1016/j.ijsolstr.2009.05.008},
+ issn = {00207683},
+ journal = {International Journal of Solids and Structures},
+ month = {September},
+ number = {18-19},
+ pages = {3357--3364},
+ title = {Modeling of high-density compaction of granular materials by the Discrete Element Method},
+ volume = {46},
+ year = {2009}
+}
+
+@article{Hartong2012a,
+ author = {Harthong, B. and Scholtès, L. and Donzé, F.-V.},
+ title = {Strength characterization of rock masses, using a coupled DEM–DFN model},
+ journal = {Geophysical Journal International},
+ volume = {191},
+ number = {2},
+ pages = {467--480},
+ year = {2012},
+ doi = {10.1111/j.1365-246X.2012.05642.x},
+ url = {http://dx.doi.org/10.1111/j.1365-246X.2012.05642.x}
+}
+
+@Article{Harthong2012b,
+ title = {Contact impingement in packings of elastic–plastic spheres, application to powder compaction},
+ journal = {International Journal of Mechanical Sciences},
+ volume = {61},
+ number = {1},
+ pages = {32–43},
+ year = {2012},
+ author = {Harthong, B. and Jerier, J.-F. and Richefeu, V. and Chareyre, B. and Doremus, P. and Imbault, D. and Donzé, F.V.}
+}
+
+@article{Hassan2010,
+ author = {Hassan, A. and Chareyre, B. and Darve, F. and Meyssonier, J. and Flin, F.},
+ title = {Microtomography-based Discrete Element Modelling of Creep in Snow},
+ journal = {Granular Matter},
+ year = {2010 (submitted)}
+}
+
+@article{Jerier2009,
+ author = {Jerier, J.-F. and Imbault, D.and Donzé, F.V. and Doremus, P.},
+ doi = {10.1007/s10035-008-0116-0},
+ journal = {Granular Matter},
+ year={2009},
+ volume={11},
+ title = {A geometric algorithm based on tetrahedral meshes to generate a dense polydisperse sphere packing},
+}
+
+@article{Jerier2010a,
+ author = {Jerier, J.-F. and Richefeu, V. and Imbault, D. and Donzé, F.V.},
+ doi = {10.1016/j.cma.2010.01.016},
+ issn = {00457825},
+ journal = {Computer Methods in Applied Mechanics and Engineering},
+ title = {Packing spherical discrete elements for large scale simulations},
+ year = {2010}
+}
+
+@article{Jerier2010b,
+ title = {Study of cold powder compaction by using the discrete element method},
+ journal = {Powder Technology},
+ volume = {In Press},
+ year = {2010},
+ doi = {10.1016/j.powtec.2010.08.056},
+ author = {Jerier, J.-F. and Hathong, B. and Richefeu, V. and Chareyre, B. and Imbault, D. and Donzé, F.-V. and Doremus, P.}
+}
+
+@article{Kozicki2005a,
+ title = {Discrete lattice model used to describe the fracture process of concrete},
+ author = {Kozicki, J.},
+ journal = {Discrete Element Group for Risk Mitigation Annual Report 1, Grenoble University of Joseph Fourier, France},
+ pages = {95--101},
+ year = {2005},
+ url={https://yade-dem.org/w/images/f/f0/Discrete_lattice_model_DEM_risk_mitigation_kozicki.pdf}
+}
+
+@article{Kozicki2006a,
+ title = {2{D} Lattice Model for Fracture in Brittle Materials},
+ author = {Kozicki, J. and Tejchman, J.},
+ journal = {Archives of Hydro-Engineering and Environmental Mechanics},
+ pages = {71--88},
+ volume = {53},
+ number = {2},
year = {2006},
- journal = {Annual Report of Discrete Element Group for Hazard Mitigation},
- url={https://yade-dem.org/w/images/a/a6/Smilauer-the_splendors_and_miseries_of_yade_design-2007.pdf}
+ url={https://yade-dem.org/w/images/5/54/Ahem_2006_kozicki.pdf}
+}
+
+@article{Kozicki2007a,
+ title = {Effect of aggregate structure on fracture process in concrete using 2D lattice model”},
+ author = {Kozicki, J. and Tejchman, J.},
+ journal = {Archives of Mechanics},
+ pages = {365--384},
+ volume = {59},
+ number = {4--5},
+ year = {2007},
+ url={https://yade-dem.org/w/images/0/09/Ams_2007_kozicki_tejchman.pdf}
}
@article{Kozicki2008,
- author = {J. Kozicki and F.V. Donzé},
+ author = {Kozicki, J. and Donzé, F.V.},
title = {A new open-source software developed for numerical simulations using discrete modeling methods},
journal = {Computer Methods in Applied Mechanics and Engineering},
volume = {197},
@@ -21,7 +318,7 @@
@article{Kozicki2009,
title={YADE-OPEN DEM: an open-source software using a discrete element method to simulate granular material},
- author={J. Kozicki and F.V. Donzé},
+ author={Kozicki, J. and Donzé, F.V.},
journal={Engineering Computations},
year={2009},
volume={26},
@@ -31,324 +328,208 @@
doi={10.1108/02644400910985170},
}
-@article{Jerier2010,
- author = {Jerier, Jean-François and Richefeu, Vincent and Imbault, Didier and Donzé, Fréderic-Victor},
- doi = {10.1016/j.cma.2010.01.016},
- issn = {00457825},
- journal = {Computer Methods in Applied Mechanics and Engineering},
- title = {Packing spherical discrete elements for large scale simulations},
- year = {2010}
-}
-
-@article{Jerier2009a,
- author = {Jerier, Jean-François and Imbault, Didier and Donzé, Fréederic-Victor and Doremus, Pierre},
- doi = {10.1007/s10035-008-0116-0},
- journal = {Granular Matter},
- year={2009},
- volume={11},
- title = {A geometric algorithm based on tetrahedral meshes to generate a dense polydisperse sphere packing},
-}
-
-@article{Harthong2009,
- author = {Harthong, B. and Jerier, J. F. and Dorémus, P. and Imbault, D. and Donzé, F. V.},
- doi = {10.1016/j.ijsolstr.2009.05.008},
- issn = {00207683},
+@Article{Lomine2013,
+ author = {Lominé, F. and Scholtès, L. and Sibille, L. and Poullain, P.},
+ title = {Modelling of fluid-solid interaction in granular media with coupled LB/DE methods: application to piping erosion},
+ journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
+ volume = {37},
+ number = {6},
+ pages = {577-596},
+ year = {2013},
+ doi = {10.1002/nag.1109}
+}
+
+@Article{Nicot2011,
+ author = {Nicot, F. and Hadda, N. and Bourrier, F. and Sibille, L. and Darve, F.},
+ title = {Failure mechanisms in granular media: a discrete element analysis},
+ journal = {Granular Matter},
+ volume = {13},
+ number = {3},
+ pages = {255-260},
+ year = {2011},
+ doi = {10.1007/s10035-010-0242-3}
+}
+
+@Article{Nicot2012,
+ author = {Nicot, F. and Sibille, L. and Darve, F.},
+ title = {Failure in rate-independent granular materials as a bifurcation toward a dynamic regime},
+ year = {2012},
+ journal = {International Journal of Plasticity},
+ volume = {29},
+ pages = {136-154},
+ doi = {10.1016/j.ijplas.2011.08.002}
+}
+
+@article{Nicot2013a,
+ title={Second-order work analysis for granular materials using a multiscale approach},
+ author={Nicot, F. and Hadda, N. and Darve, F.},
+ journal={International Journal for Numerical and Analytical Methods in Geomechanics},
+ year={2013},
+ doi={10.1002/nag.2175}
+}
+
+@article{Nicot2013b,
+ title = {On the definition of the stress tensor in granular media},
journal = {International Journal of Solids and Structures},
- month = {September},
- number = {18-19},
- pages = {3357--3364},
- title = {Modeling of high-density compaction of granular materials by the Discrete Element Method},
- volume = {46},
- year = {2009}
+ year = {2013},
+ doi = {10.1016/j.ijsolstr.2013.04.001},
+ url = {http://www.sciencedirect.com/science/article/pii/S0020768313001492},
+ author = {Nicot, F. and Hadda, N. and Guessasma, M. and Fortin, J. and Millet, O.}
+}
+
+@article{Puckett2011,
+ title = {Local origins of volume fraction fluctuations in dense granular materials},
+ author = {Puckett, J.G. and Lechenault, F. and Daniels, K.E.},
+ journal = {Physical Review E},
+ volume = {83},
+ issue = {4},
+ pages = {041301},
+ year = {2011},
+ doi = {10.1103/PhysRevE.83.041301},
+ url = {http://link.aps.org/doi/10.1103/PhysRevE.83.041301}
+}
+
+@article{Sayeed2011,
+ title={Strength and Deformation Characteristics of Granular Materials under Extremely Low to High Confining Pressures in Triaxial Compression},
+ author={Sayeed, M.A. and Suzuki, K. and Rahman, M.M. and Mohamad, W.H.W. and Razlan, M.A. and Ahmad, Z. and Thumrongvut, J. and Seangatith, S. and Sobhan, MA and Mofiz, SA and others},
+ volume={11},
+ number={4},
+ journal={International Journal of Civil & Environmental Engineering IJCEE-IJENS},
+ year={2011}
}
@article{Scholtes2009a,
author = {Scholtès, L. and Chareyre, B. and Nicot, F. and Darve, F.},
- doi = {10.1016/j.ijengsci.2008.07.002},
- issn = {},
+ title = {Micromechanics of granular materials with capillary effects},
journal = {International Journal of Engineering Science},
- month = {},
+ volume = {47},
number = {1},
pages = {64--75},
- title = {Micromechanics of granular materials with capillary effects},
- volume = {47},
- year = {2009}
+ year = {2009},
+ doi = {10.1016/j.ijengsci.2008.07.002},
+ url = {http://dx.doi.org/10.1016/j.ijengsci.2008.07.002}
}
@article{Scholtes2009b,
author = {Scholtès, L. and Hicher, P.-Y. and Chareyre, B. and Nicot, F. and Darve, F.},
- doi = {10.1002/nag.767},
- issn = {},
+ title = {On the capillary stress tensor in wet granular materials},
journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
- month = {},
+ volume = {33},
number = {10},
pages = {1289--1313},
- title = {On the capillary stress tensor in wet granular materials},
- volume = {33},
- year = {2009}
+ year = {2009},
+ doi = {10.1002/nag.767},
+ url = {http://arxiv.org/abs/1105.1013}
}
@article{Scholtes2009c,
author = {Scholtès, L. and Chareyre, B. and Nicot, F. and Darve, F.},
- doi = {},
- issn = {1526-1506},
+ title = {Discrete modelling of capillary mechanisms in multi-phase granular media},
journal = {Computer Modeling in Engineering and Sciences},
- month = {},
+ volume = {52},
number = {3},
pages = {297--318},
- title = {Discrete modelling of capillary mechanisms in multi-phase granular media},
+ year = {2009},
+ url = {http://arxiv.org/abs/1203.1234}
+}
+
+@Article{Scholtes2010,
+ author = {Scholtès, L. and Hicher, P.-Y. and Sibille, L.},
+ title = {Multiscale approaches to describe mechanical responses induced by particle removal in granular materials},
+ journal = {Comptes Rendus Mécanique},
+ volume = {338},
+ number = {10-11},
+ pages = {627--638},
+ year = {2010},
+ doi = {10.1016/j.crme.2010.10.003},
+ url = {http://dx.doi.org/10.1016/j.crme.2010.10.003}
+}
+
+@article{Scholtes2011,
+ author={Scholtès, L. and Donzé, F.V. and Khanal, M.},
+ title={Scale effects on strength of geomaterials, case study: Coal},
+ journal={Journal of the Mechanics and Physics of Solids},
+ volume={59},
+ number={5},
+ pages={1131--1146},
+ year={2011},
+ doi={10.1016/j.jmps.2011.01.009},
+ url = {http://dx.doi.org/10.1016/j.jmps.2011.01.009}
+}
+
+@article{Scholtes2012,
+ author = {Scholtès, L. and Donzé, F.V.},
+ title = {Modelling progressive failure in fractured rock masses using a 3D Discrete Element Method},
+ journal = {International Journal of Rock Mechanics and Mining Sciences},
volume = {52},
- year = {2009}
-}
-
-@article{Chen2007,
- title = {Prediction/Verification of Particle Motion in One Dimension with the Discrete-Element Method},
- author = {Feng Chen and Eric. C. Drumm and Georges Guiochon},
- doi = {10.1061/(ASCE)1532-3641(2007)7:5(344)},
- journal = {International Journal of Geomechanics, ASCE},
- pages = {344--352},
- volume = {7},
- number = {5},
- year = {2007}
-}
-
-@article{Kozicki2007a,
- title = {Effect of aggregate structure on fracture process in concrete using 2D lattice model”},
- author = {J. Kozicki and J. Tejchman},
- journal = {Archives of Mechanics},
- pages = {365--384},
- volume = {59},
- number = {4--5},
- year = {2007},
- url={https://yade-dem.org/w/images/0/09/Ams_2007_kozicki_tejchman.pdf}
-}
-
-@article{Kozicki2006a,
- title = {2D Lattice Model for Fracture in Brittle Materials},
- author = {J. Kozicki and J. Tejchman},
- journal = {Archives of Hydro-Engineering and Environmental Mechanics},
- pages = {71--88},
- volume = {53},
- number = {2},
- year = {2006},
- url={https://yade-dem.org/w/images/5/54/Ahem_2006_kozicki.pdf}
-}
-
-@article{Kozicki2005a,
- title = {Discrete lattice model used to describe the fracture process of concrete},
- author = {J. Kozicki},
- journal = {Discrete Element Group for Risk Mitigation Annual Report 1, Grenoble University of Joseph Fourier, France},
- pages = {95--101},
- year = {2005},
- url={https://yade-dem.org/w/images/f/f0/Discrete_lattice_model_DEM_risk_mitigation_kozicki.pdf}
-}
-
-@article{ Hassan2010,
- author = {A. Hassan and B. Chareyre and F. Darve and J. Meyssonier and F. Flin},
- title = {Microtomography-based Discrete Element Modelling of Creep in Snow},
- journal = {Granular Matter},
- year = {2010 (submitted)}
-}
-
-@article{Duriez2011,
- author = {J. Duriez and F.Darve and F.-V.Donze},
- title = {A discrete modeling-based constitutive relation for infilled rock joints},
- year = {2011},
- journal = {International Journal of Rock Mechanics & Mining Sciences},
- volume = {48},
- number = {3},
- pages = {458--468},
- doi = {10.1016/j.ijrmms.2010.09.008}
-}
-
-@article{Dang2010a,
- author = {H. K. Dang and M. A. Meguid},
- title = {Algorithm to Generate a Discrete Element Specimen with Predefined Properties},
- publisher = {ASCE},
- year = {2010},
- journal = {International Journal of Geomechanics},
- volume = {10},
- number = {2},
- pages = {85-91},
- keywords = {Discrete elements; Geotechnical engineering; Grain size; Porosity; Algorithms; Granular materials},
- doi = {10.1061/(ASCE)GM.1943-5622.0000028}
-}
-
-@article{ Jerier2010b,
- title = {Study of cold powder compaction by using the discrete element method},
- journal = {Powder Technology},
- volume = {In Press},
- year = {2010},
- doi = {10.1016/j.powtec.2010.08.056},
- author = {J.-F. Jerier and B. Hathong and V. Richefeu and B. Chareyre and D. Imbault and F.-V. Donze and P. Doremus}
-}
-
-@article{Dang2010b,
- title = {Evaluating the performance of an explicit dynamic relaxation technique in analyzing non-linear geotechnical engineering problems},
- journal = {Computers and Geotechnics},
- volume = {37},
- number = {1-2},
- pages = {125 - 131},
- year = {2010},
- doi = {10.1016/j.compgeo.2009.08.004},
- author = {Hoang K. Dang and Mohamed A. Meguid}
-}
-
-@article{Chen2011a,
- title={Coupled discrete element and finite volume solution of two classical soil mechanics problems},
- author={Chen, F. and Drumm, E. and Guiochon G.},
- journal={Computers and Geotechnics},
- year={2011},
- doi = {10.1016/j.compgeo.2011.03.009},
- url={http://www.sciencedirect.com/science/article/pii/S0266352X11000504},
- publisher={Elsevier}
+ pages = {18--30},
+ year = {2012},
+ doi = {10.1016/j.ijrmms.2012.02.009},
+ url = {http://dx.doi.org/10.1016/j.ijrmms.2012.02.009}
+}
+
+@article{Scholtes2013,
+ author = {Scholtès, L. and Donzé, F.V.},
+ title = {A {DEM} model for soft and hard rocks: Role of grain interlocking on strength},
+ journal = {Journal of the Mechanics and Physics of Solids},
+ volume = {61},
+ number = {2},
+ pages = {352--369},
+ year = {2013},
+ doi = {10.1016/j.jmps.2012.10.005},
+ url = {http://dx.doi.org/10.1016/j.jmps.2012.10.005}
}
@article{Shiu2008,
- title={Compaction process in concrete during missile impact: a DEM analysis},
- author={Shiu, W. and Donzé, FV and Daudeville, L.},
- journal={Computers and Concrete},
- volume={5},
- number={4},
- pages={329--342},
- year={2008}
+ title={Compaction process in concrete during missile impact: a DEM analysis},
+ author={Shiu, W. and Donzé, F.V. and Daudeville, L.},
+ journal={Computers and Concrete},
+ volume={5},
+ number={4},
+ pages={329--342},
+ year={2008}
}
@article{Shiu2009,
- title={Discrete element modelling of missile impacts on a reinforced concrete target},
- author={Shiu, W. and Donze, F.V. and Daudeville, L.},
- journal={International Journal of Computer Applications in Technology},
- volume={34},
- number={1},
- pages={33--41},
- year={2009},
-}
-
-@article{Tran2011a,
- title={A discrete element model of concrete under high triaxial loading},
- author={Tran, VT and Donzé, F.V. and Marin, P.},
- journal={Cement and Concrete Composites},
- year={2011},
- publisher={Elsevier}
-}
-
-@article{Sayeed2011,
- title={Strength and Deformation Characteristics of Granular Materials under Extremely Low to High Confining Pressures in Triaxial Compression},
- author={Sayeed, M.A. and Suzuki, K. and Rahman, M.M. and Mohamad, W.H.W. and Razlan, M.A. and Ahmad, Z. and Thumrongvut, J. and Seangatith, S. and Sobhan, MA and Mofiz, SA and others},
- volume={11},
- number={4},
- journal={International Journal of Civil & Environmental Engineering IJCEE-IJENS},
- year={2011}
-}
-
-@article{Duriez2011,
- title={Incrementally non-linear plasticity applied to rock joint modelling},
- author={Duriez, J. and Darve, F. and Donzé, F.V.},
- journal={International Journal for Numerical and Analytical Methods in Geomechanics},
- year={2011},
- doi = {10.1002/nag.1105}
-}
-
-@article{Donze2008,
- title={Impacts on cohesive frictional geomaterials},
- author={Donzé, F.V.},
- journal={European Journal of Environmental and Civil Engineering},
- volume={12},
- number={7-8},
- pages={967--985},
- year={2008}
-}
-
-@article{Scholtes2011,
- title={Scale effects on strength of geomaterials, case study: Coal},
- author={Scholtès, L. and Donzé, F.V. and Khanal, M.},
- journal={Journal of the Mechanics and Physics of Solids},
- year={2011}
+ title={Discrete element modelling of missile impacts on a reinforced concrete target},
+ author={Shiu, W. and Donzé, F.V. and Daudeville, L.},
+ journal={International Journal of Computer Applications in Technology},
+ volume={34},
+ number={1},
+ pages={33--41},
+ year={2009},
+}
+
+@article{Smilauer2006,
+ title = {The splendors and miseries of Yade design},
+ author = {Václav Šmilauer},
+ year = {2006},
+ journal = {Annual Report of Discrete Element Group for Hazard Mitigation},
+ url={https://yade-dem.org/w/images/a/a6/Smilauer-the_splendors_and_miseries_of_yade_design-2007.pdf}
}
@article{Tejchman2011,
- title={Comparative Modelling of Shear Zone Patterns in Granular Bodies with Finite and Discrete Element Model},
- author={Tejchman, J.},
- journal={Advances in Bifurcation and Degradation in Geomaterials},
- pages={255--260},
- year={2011}
-}
-
-@article{Favier2009a,
- title={Predicting the drag coefficient of a granular flow using the discrete element method},
- author={Favier, L. and Daudon, D. and Donzé, F.V. and Mazars, J.},
- journal={Journal of Statistical Mechanics: Theory and Experiment},
- volume={2009},
- pages={P06012},
- year={2009}
-}
-
-@article {Chareyre2012a,
- author = {Chareyre, B. and Cortis, A. and Catalano, E. and Barthélemy, E.},
- title = {Pore-Scale Modeling of Viscous Flow and Induced Forces in Dense Sphere Packings},
- journal = {Transport in Porous Media},
- publisher = {Springer Netherlands},
- volume={92},
- number={2},
- year={2012},
- pages = {473-493},
- url = {http://dx.doi.org/10.1007/s11242-011-9915-6},
- doi = {10.1007/s11242-011-9915-6}
-}
-
-@Article{ Duriez2011,
- author = {J. Duriez and F. Darve and F.-V. Donze},
- title = {A discrete modeling-based constitutive relation for infilled rock joints},
- year = {2011},
- journal = {International Journal of Rock Mechanics \& Mining Sciences},
- volume = {48},
- number = {3},
- pages = {458--468},
- doi = {10.1016/j.ijrmms.2010.09.008}
-}
-
-@Article{Nicot2011a,
- author = {Nicot, F. and Hadda, N. and Bourrier, F. and Sibille, L. and Darve, F.},
- title = {Failure mechanisms in granular media: a discrete element analysis},
- year = {2011},
- journal = {Granular Matter},
- volume = {13},
- number = {3},
- pages = {255-260},
- doi = {10.1007/s10035-010-0242-3}
-}
-
-@Article{Nicot2012,
- author = {Nicot, F. and Sibille, L. and Darve, F.},
- title = {Failure in rate-independent granular materials as a bifurcation toward a dynamic regime},
- year = {2012},
- journal = {International Journal of Plasticity},
- volume = {29},
- pages = {136-154},
- doi = {10.1016/j.ijplas.2011.08.002}
-}
-
-@Article{Scholtes2010,
- author = {Scholtès, L. and Hicher, P.Y. and Sibille, L.},
- title = {Multiscale approaches to describe mechanical responses induced by particle removal in granular materials},
- year = {2010},
- journal = {Comptes Rendus Mécanique (CRAS)},
- volume = {338},
- number = {10-11},
- pages = {627-638},
- doi = {10.1016/j.crme.2010.10.003}
-}
-
-@Article{Lomine2012,
- author = {Lominé, F. and Scholtès, L. and Sibille, L. and Poullain, P.},
- title = {Modelling of fluid-solid interaction in granular media with coupled LB/DE methods: application to piping erosion},
- year = {2012},
- journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
- doi = {10.1002/nag.1109}
+ title={Comparative Modelling of Shear Zone Patterns in Granular Bodies with Finite and Discrete Element Model},
+ author={Tejchman, J.},
+ journal={Advances in Bifurcation and Degradation in Geomaterials},
+ pages={255--260},
+ year={2011}
+}
+
+@article{Thoeni2013,
+ author = "K. Thoeni and C. Lambert and A. Giacomini and S.W. Sloan",
+ doi = "10.1016/j.compgeo.2012.10.014",
+ issn = "0266-352X",
+ journal = "Computers and Geotechnics",
+ pages = "158--169",
+ title = "{Discrete modelling of hexagonal wire meshes with a stochastically distorted contact model}",
+ url = "http://www.sciencedirect.com/science/article/pii/S0266352X12002121";,
+ volume = "49",
+ year = "2013"
}
@article{Tong2012,
- author = {{Tong, A.-T.} and {Catalano, E.} and {Chareyre, B.}},
+ author = {Tong, A.-T. and Catalano, E. and Chareyre, B.},
title = {Pore-Scale Flow Simulations: Model Predictions Compared with Experiments on Bi-Dispersed Granular Assemblies},
doi= {10.2516/ogst/2012032},
url= {http://dx.doi.org/10.2516/ogst/2012032},
@@ -356,85 +537,25 @@
year = {2012}
}
-@article{Tran2012a,
- title={An Algorithm for the Propagation of Uncertainty in Soils using the Discrete Element Method},
- author={Tran, V.D.H. and Meguid, M.A. and Chouinard, L.E.},
- journal = {The Electronic Journal of Geotechnical Engineering},
- url={http://www.ejge.com/2012/Ppr12.283alr.pdf},
- year= {2012}
-}
-
-@article{Boon2012,
- title = {A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method},
- journal = {Computers and Geotechnics},
- volume = {44},
- number = {0},
- pages = {73 - 82},
- year = {2012},
- doi = {10.1016/j.compgeo.2012.03.012},
- url = {http://www.sciencedirect.com/science/article/pii/S0266352X12000535},
- author = {C.W. Boon and G.T. Houlsby and S. Utili}
-}
-
-@article {Hartong2012a,
- author = {Harthong, Barthélémy and Scholtès, Luc and Donzé, Frédéric-Victor},
- title = {Strength characterization of rock masses, using a coupled DEM–DFN model},
- journal = {Geophysical Journal International},
- volume = {191},
- number = {2},
- url = {http://dx.doi.org/10.1111/j.1365-246X.2012.05642.x},
- doi = {10.1111/j.1365-246X.2012.05642.x},
- pages = {467--480},
- year = {2012}
-}
-
-@Article{ Harthong2012b,
- title = {Contact impingement in packings of elastic–plastic spheres, application to powder compaction},
- journal = {International Journal of Mechanical Sciences},
- volume = {61},
- number = {1},
- pages = {32–43},
- year = {2012},
- author = {Barthélémy Harthong and Jean-François Jérier and Vincent Richefeu and Bruno Chareyre and Pierre Dorémus and Didier Imbault and Frédéric-Victor Donzé}
-}
-
-@article{Puckett2011,
- title = {Local origins of volume fraction fluctuations in dense granular materials},
- author = {Puckett, James G. and Lechenault, Frédéric and Daniels, Karen E.},
- journal = {Phys. Rev. E},
- volume = {83},
- issue = {4},
- pages = {041301},
- numpages = {8},
- year = {2011},
- month = {Apr},
- doi = {10.1103/PhysRevE.83.041301},
- url = {http://link.aps.org/doi/10.1103/PhysRevE.83.041301},
-}
-
-@article{Chen2012,
- author = {Chen, Jingsong and Huang, Baoshan and Chen, Feng and Shu, Xiang},
- title = {Application of discrete element method to Superpave gyratory compaction},
- journal = {Road Materials and Pavement Design},
- volume = {13},
- number = {3},
- pages = {480-500},
- year = {2012},
- doi = {10.1080/14680629.2012.694160},
- URL = {http://www.tandfonline.com/doi/abs/10.1080/14680629.2012.694160}
-}
-
-@article{Nicot2013,
- title={Second-order work analysis for granular materials using a multiscale approach},
- author={Nicot, François and Hadda, Nejib and Darve, Félix},
- journal={International Journal for Numerical and Analytical Methods in Geomechanics},
- year={2013},
- doi={10.1002/nag.2175}
+@article{Tran2011,
+ title={A discrete element model of concrete under high triaxial loading},
+ author={Tran, V.T. and Donzé, F.V. and Marin, P.},
+ journal={Cement and Concrete Composites},
+ year={2011},
+ publisher={Elsevier}
+}
+
+@article{Tran2012,
+ title={An Algorithm for the Propagation of Uncertainty in Soils using the Discrete Element Method},
+ author={Tran, V.D.H. and Meguid, M.A. and Chouinard, L.E.},
+ journal = {The Electronic Journal of Geotechnical Engineering},
+ url={http://www.ejge.com/2012/Ppr12.283alr.pdf},
+ year= {2012}
}
@article{Tran2013,
title={A finite--discrete element framework for the 3D modeling of geogrid--soil interaction under pullout loading conditions},
- author={Tran, VDH and Meguid, MA and Chouinard, LE},
+ author={Tran, V.D.H. and Meguid, M.A. and Chouinard, L.E.},
journal={Geotextiles and Geomembranes},
volume={37},
pages={1-9},
@@ -445,125 +566,44 @@
@article{Tran2012c,
title={Discrete Element and Experimental Investigations of the Earth Pressure Distribution on Cylindrical Shafts},
- author={Tran, Viet DH and Meguid, Mohamed A and Chouinard, Luc E},
+ author={Tran, V.D.H. and Meguid, M.A. and Chouinard, L.E.},
journal={International Journal of Geomechanics},
year={2012},
publisher={American Society of Civil Engineers},
doi = {10.1061/(ASCE)GM.1943-5622.0000277}
}
-@article{Epif2012a,
- title={Proizvodstvo kuskovogo torfa, ekstrudirovanie, forma zakhodnoi i kalibriruyushchei chasti fil'ery matritsy, metod diskretnykh elementov [RUS]},
- author={Epifantsev, K. and Mikhailov, A. and Gladky, A.},
- journal={Mining informational and analytical bulletin (scientific and technical journal)},
- year={2012},
- publisher={The joint-stock company Gornaya Kniga},
- number = {3},
- pages = {212-219},
- issn = {0236-1493},
-}
-
-
-@article {Duriez2013,
- author = {Duriez, Jérôme and Darve, Félix and Donzé, Frédéric-Victor},
- title = {Incrementally non-linear plasticity applied to rock joint modelling},
- journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
- volume = {37},
- number = {5},
- url = {http://dx.doi.org/10.1002/nag.1105},
- doi = {10.1002/nag.1105},
- pages = {453--477},
- year = {2013}
-}
-
-@article{Nicot2013,
- title = {On the definition of the stress tensor in granular media},
- journal = {International Journal of Solids and Structures},
- year = {2013},
- doi = {10.1016/j.ijsolstr.2013.04.001},
- url = {http://www.sciencedirect.com/science/article/pii/S0020768313001492},
- author = {François Nicot and Nejib Hadda and Mohamed Guessasma and Jerome Fortin and Olivier Millet}
-}
-
-@article{Gusenbauer2012,
- title={Self-organizing magnetic beads for biomedical applications},
- author={Gusenbauer, Markus and Kovacs, Alexander and Reichel, Franz and Exl, Lukas and Bance, Simon and {\"O}zelt, Harald and Schrefl, Thomas},
- journal={Journal of Magnetism and Magnetic Materials},
- volume={324},
- number={6},
- pages={977--982},
- year={2012},
- publisher={Elsevier}
-}
-
@article{Bourrier2013,
title = {Discrete modeling of granular soils reinforcement by plant roots},
journal = {Ecological Engineering},
- author={Bourrier, F and Kneib, F and Chareyre, B and Fourcaud, T},
+ author={Bourrier, F. and Kneib, F. and Chareyre, B. and Fourcaud, T.},
year = {2013},
doi = {10.1016/j.ecoleng.2013.05.002},
url = {http://dx.doi.org/10.1016/j.ecoleng.2013.05.002}
}
-@article{Grujicic2013,
- title={Discrete Element Modeling and Analysis of Structural Collapse/Survivability of a Building Subjected to Improvised Explosive Device (IED) Attack},
- author={Grujicic, M and Snipes, JS and Ramaswami, S and Yavari, R},
- journal={Advances in Materials Science and Applications},
- year={2013},
- volume={2},
- number={1},
- pages={9--24}
-}
-
-@article{Scholtès2013,
- title = {A \{DEM\} model for soft and hard rocks: Role of grain interlocking on strength},
- journal = {Journal of the Mechanics and Physics of Solids},
- volume = {61},
- number = {2},
- pages = {352--369},
- year = {2013},
- doi = {10.1016/j.jmps.2012.10.005},
- url = {http://www.sciencedirect.com/science/article/pii/S0022509612002268},
- author = {Luc Scholtès and Frédéric-Victor Donzé}
-}
-
-@article{Boon2012,
- title = {A new contact detection algorithm for three-dimensional non-spherical particles},
- journal = {Powder Technology},
- year = {2013},
- doi = {10.1016/j.powtec.2012.12.040},
- url = {http://www.sciencedirect.com/science/article/pii/S003259101200839X},
- author = {C.W. Boon and G.T. Houlsby and S. Utili}
-}
-
-@article{Favier2012,
- title={Rigid obstacle impacted by a supercritical cohesive granular flow using a 3D discrete element model},
- author={Favier, L. and Daudon, D. and Donzé, F.V.},
- journal={Cold Regions Science and Technology},
- year={2012},
- volume = {85},
- pages = {232--241},
- url = {http://dx.doi.org/10.1016/j.coldregions.2012.09.010}
-}
-
-@Article{Catalano2013a,
- title = {Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects},
- author = {E. Catalano and B. Chareyre and E. Barthélémy},
- journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
- year = {2013},
- doi = {10.1002/nag.2198},
- note = {in press, arxiv version available},
- url = {http://arxiv.org/pdf/1304.4895.pdf}
-}
-
-@article{Thoeni2013,
- author = "K. Thoeni and C. Lambert and A. Giacomini and S.W. Sloan",
- doi = "10.1016/j.compgeo.2012.10.014",
- issn = "0266-352X",
- journal = "Computers and Geotechnics",
- pages = "158--169",
- title = "{Discrete modelling of hexagonal wire meshes with a stochastically distorted contact model}",
- url = "http://www.sciencedirect.com/science/article/pii/S0266352X12002121";,
- volume = "49",
- year = "2013"
-}
\ No newline at end of file
+@article{Epifancev2013,
+ author={Epifancev, K. and Nikulin, A. and Kovshov, S. and Mozer, S. and Brigadnov, I.},
+ title={Modeling of Peat Mass Process Formation Based on 3D Analysis of the Screw Machine by the Code YADE},
+ journal={American Journal of Mechanical Engineering},
+ volume={1},
+ number={3},
+ pages={73--75},
+ year={2013},
+ url={http://pubs.sciepub.com/ajme/1/3/3},
+ doi={10.12691/ajme-1-3-3}
+}
+
+@article{Hadda2013,
+ year={2013},
+ journal={Granular Matter},
+ volume={15},
+ number={2},
+ doi={10.1007/s10035-013-0402-3},
+ title={Micromechanical analysis of second order work in granular media},
+ url={http://dx.doi.org/10.1007/s10035-013-0402-3},
+ author={Hadda, Nejib and Nicot, François and Bourrier, Franck and Sibille, Luc and Radjai, Farhang and Darve, Félix},
+ pages={221--235}
+}
+
+
=== modified file 'doc/yade-conferences.bib'
--- doc/yade-conferences.bib 2013-03-05 10:48:14 +0000
+++ doc/yade-conferences.bib 2013-07-29 17:38:12 +0000
@@ -505,3 +505,88 @@
year={2012},
month={October}
}
+
+
+@InProceedings{ Harthong2013,
+ title = {Directional plastic flow and fabric dependencies in granular materials},
+ author = {Barthélémy Harthong and Richard G Wan},
+ Booktitle = {{Powders and Grains 2013: Proceedings of the 6th International Conference on Micromechanics of Granular Media. AIP Conference Proceedings}},
+ volume = {1542},
+ year = {2009},
+ number = {1},
+ month = {July},
+ doi = {http://dx.doi.org/10.1063/1.4811900},
+ pages = {193–197},
+ address = {Sydney, Australia},
+ url = {https://www.yade-dem.org/publi/APC000193.pdf}
+}
+
+@InProceedings{ Hilton2013,
+ author = {J. E. Hilton and P. W. Cleary and A. Tordesillas},
+ title = {Unitary stick-slip motion in granular beds},
+ publisher = {AIP},
+ year = {2013},
+ Booktitle = {{Powders and Grains 2013: Proceedings of the 6th International Conference on Micromechanics of Granular Media. AIP Conference Proceedings}},
+ volume = {1542},
+ number = {1},
+ pages = {843–846},
+ address = {Sydney, Australia},
+ url = {http://scitation.aip.org/getpdf/servlet/GetPDFServlet?filetype=pdf&id=APCPCS001542000001000843000001&idtype=cvips&doi=10.1063/1.4812063&prog=normal},
+ doi = {10.1063/1.4812063}
+}
+
+@InProceedings{ Catalano2013b,
+ author = {Catalano E. and Chareyre B. and Barthélémy E.},
+ title = {DEM-PFV analysis of solid-fluid transition in granular sediments under the action of waves},
+ publisher = {AIP},
+ year = {2013},
+ Booktitle = {{Powders and Grains 2013: Proceedings of the 6th International Conference on Micromechanics of Granular Media. AIP Conference Proceedings}},
+ volume = {1542},
+ number = {1},
+ pages = {1063–1066},
+ address = {Sydney, Australia},
+ url = {http://link.aip.org/link/?APC/1542/1063/1},
+ doi = {10.1063/1.4812118}
+}
+
+@InProceedings{ Guo2013,
+ author = {Ning Guo and Jidong Zhao},
+ title = {A hierarchical model for cross-scale simulation of granular media},
+ publisher = {AIP},
+ year = {2013},
+ Booktitle = {{Powders and Grains 2013: Proceedings of the 6th International Conference on Micromechanics of Granular Media. AIP Conference Proceedings}},
+ volume = {1542},
+ number = {1},
+ pages = {1222–1225},
+ address = {Sydney, Australia},
+ url = {https://www.yade-dem.org/publi/APC001222.pdf},
+ doi = {10.1063/1.4812158}
+}
+
+@InProceedings{ Hadda2013,
+ author = {Nejib Hadda and François Nicot and Luc Sibille and Farhang Radjai and Antoinette Tordesillas and Félix Darve},
+ title = {A multiscale description of failure in granular materials},
+ publisher = {AIP},
+ year = {2013},
+ Booktitle = {{Powders and Grains 2013: Proceedings of the 6th International Conference on Micromechanics of Granular Media. AIP Conference Proceedings}},
+ volume = {1542},
+ number = {1},
+ pages = {585–588},
+ address = {Sydney, Australia},
+ url = {https://www.yade-dem.org/publi/APC000585.pdf},
+ doi = {10.1063/1.4811999}
+}
+
+@InProceedings{ Kozicki2013,
+ author = {Jan Kozicki and Jacek Tejchman and Danuta Lesniewska},
+ title = {Study of some micro-structural phenomena in granular shear zones},
+ publisher = {AIP},
+ year = {2013},
+ Booktitle = {{Powders and Grains 2013: Proceedings of the 6th International Conference on Micromechanics of Granular Media. AIP Conference Proceedings}},
+ volume = {1542},
+ number = {1},
+ pages = {495–498},
+ address = {Sydney, Australia},
+ url = {https://www.yade-dem.org/publi/APC000495.pdf},
+ doi = {10.1063/1.4811976}
+}
=== added directory 'examples/FluidCouplingPFV'
=== added file 'examples/FluidCouplingPFV/oedometer.py'
--- examples/FluidCouplingPFV/oedometer.py 1970-01-01 00:00:00 +0000
+++ examples/FluidCouplingPFV/oedometer.py 2013-07-29 17:36:38 +0000
@@ -0,0 +1,165 @@
+# -*- coding: utf-8 -*-
+#*************************************************************************
+# Copyright (C) 2010 by Bruno Chareyre *
+# bruno.chareyre_at_grenoble-inp.fr *
+# *
+# This program is free software; it is licensed under the terms of the *
+# GNU General Public License v2 or later. See file LICENSE for details. *
+#*************************************************************************/
+
+## Example script for using the DEM-PFV coupling introduced with E. Catalano, as reported in:
+## * [Chareyre2012a] Chareyre, B., Cortis, A., Catalano, E., Barthélemy, E. (2012), Pore-scale modeling of viscous flow and induced forces in dense sphere packings. Transport in Porous Media (92), pages 473-493. DOI 10.1007/s11242-011-9915-6
+## http://dx.doi.org/10.1007/s11242-011-9915-6
+## * [Catalano2013a] Catalano, E., Chareyre, B., Barthélémy, E. (2013), Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects. International Journal for Numerical and Analytical Methods in Geomechanics. DOI 10.1002/nag.2198
+## http://arxiv.org/pdf/1304.4895.pdf
+## Also used in:
+## * Tong et al.2012 (http://dx.doi.org/10.2516/ogst/2012032)
+## * Sari et al 2011 (http://geo.hmg.inpg.fr/~chareyre/pubs/SariChareyreCatalanoPhilippeVincens_Particles2011.pdf)
+
+
+## The DEM-PFV is applied here to 1D consolidation (oedometer test). The example includes the determination of oedometer modulus Ee and permeability K.
+## The 1D consolidation is simulated as a coupled problem and the analytical solution corresponding to the abovementionned Ee and K is used for comparison.
+## See triax-tutorial/script-session1.py for more detailed explanations of the packing generation procedure.
+
+## ______________ First section, similar to triax-tutorial/script-session1.py _________________
+from yade import pack
+
+num_spheres=1000# number of spheres
+young=1e6
+compFricDegree = 3 # initial contact friction during the confining phase
+finalFricDegree = 30 # contact friction during the deviatoric loading
+mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing
+
+O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
+O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
+walls=aabbWalls([mn,mx],thickness=0,material='walls')
+wallIds=O.bodies.append(walls)
+
+sp=pack.SpherePack()
+sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
+sp.toSimulation(material='spheres')
+
+triax=TriaxialStressController(
+ maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
+ finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
+ thickness = 0,
+ stressMask = 7,
+ max_vel = 0.005,
+ internalCompaction=True, # If true the confining pressure is generated by growing particles
+)
+
+newton=NewtonIntegrator(damping=0.2)
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
+ [Ip2_FrictMat_FrictMat_FrictPhys()],
+ [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"
+ ),
+ FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
+ GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
+ triax,
+ newton
+]
+
+triax.goal1=triax.goal2=triax.goal3=10000
+
+while 1:
+ O.run(1000, True)
+ unb=unbalancedForce()
+ if unb<0.001 and abs(10000-triax.meanStress)/10000<0.001:
+ break
+
+setContactFriction(radians(finalFricDegree))
+
+## ______________ Oedometer section _________________
+
+#A. Check bulk modulus of the dry material from load/unload cycles
+triax.stressMask=2
+triax.goal1=triax.goal3=0
+
+triax.internalCompaction=False
+triax.wall_bottom_activated=False
+#load
+triax.goal2=11000; O.run(2000,1)
+#unload
+triax.goal2=10000; O.run(2000,1)
+#load
+triax.goal2=11000; O.run(2000,1)
+e22=triax.strain[1]
+#unload
+triax.goal2=10000; O.run(2000,1)
+
+e22=e22-triax.strain[1]
+modulus = 1000./abs(e22)
+
+#B. Activate flow engine and set boundary conditions in order to get permeability
+flow.dead=0
+flow.defTolerance=0.3
+flow.meshUpdateInterval=200
+flow.useSolver=3
+flow.permeabilityFactor=1
+flow.viscosity=10
+flow.bndCondIsPressure=[0,0,1,1,0,0]
+flow.bndCondValue=[0,0,1,0,0,0]
+flow.boundaryUseMaxMin=[0,0,0,0,0,0]
+O.dt=0.1e-3
+O.dynDt=False
+
+O.run(1,1)
+Qin = flow.getBoundaryFlux(2)
+Qout = flow.getBoundaryFlux(3)
+permeability = abs(Qin)/1.e-4 #size is one, we compute K=V/∇H
+print "Qin=",Qin," Qout=",Qout," permeability=",permeability
+
+#C. now the oedometer test, drained at the top, impermeable at the bottom plate
+flow.bndCondIsPressure=[0,0,0,1,0,0]
+flow.bndCondValue=[0,0,0,0,0,0]
+newton.damping=0
+
+#we want the theoretical value from Terzaghi's solution
+#keep in mind that we are not in an homogeneous material and the smal strain
+#assumption is not verified => we don't expect perfect match
+#there can be also an overshoot of pressure in the very beginning due to dynamic effects
+Cv=permeability*modulus/1e4
+zeroTime=O.time
+zeroe22 = triax.strain[1]
+dryFraction=0.05 #the top layer is affected by drainage on a certain depth, we account for it here
+drye22 = 1000/modulus*dryFraction
+wetHeight=1*(1-dryFraction)
+
+def consolidation(Tv): #see your soil mechanics handbook...
+ U=1
+ for k in range(50):
+ M=pi/2*(2*k+1)
+ U=U-2/M**2*exp(-M**2*Tv)
+ return U
+
+triax.goal2=11000
+
+
+
+from yade import plot
+
+## a function saving variables
+def history():
+ plot.addData(e22=triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=triax.stress(3)[1]-10000)
+ #plot.addData(e22=triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5)))
+
+O.engines=O.engines+[PyRunner(iterPeriod=200,command='history()',label='recorder')]
+##make nice animations:
+#O.engines=O.engines+[PyRunner(iterPeriod=200,command='flow.saveVtk()')]
+
+from yade import plot
+plot.plots={'t':('e22','e22_theory',None,'s22','p')}
+plot.plot()
+O.saveTmp()
+O.timingEnabled=1
+from yade import timing
+print "starting oedometer simulation"
+O.run(200,1)
+timing.stats()
+
+## Make more steps to see the convergence to the stationnary solution
=== added file 'examples/LudingPM.py'
--- examples/LudingPM.py 1970-01-01 00:00:00 +0000
+++ examples/LudingPM.py 2013-07-09 12:02:21 +0000
@@ -0,0 +1,67 @@
+#!/usr/bin/env python
+# encoding: utf-8
+from yade import utils, plot
+o = Omega()
+
+fr = 0.5;rho=2000
+
+o.dt = 0.000002
+
+r1 = 10.0e-2
+r2 = 10.0e-2
+
+mat1 = O.materials.append(LudingMat(frictionAngle=fr,density=rho,k1=0.2, kp=0.9, kc=0.1,PhiF=0.01, G0 = 0.0))
+
+id1 = O.bodies.append(sphere(center=[0,0,0],radius=r1,material=mat1,fixed=True))
+id2 = O.bodies.append(sphere(center=[0,0,(r1 + r2)],radius=r2,material=mat1,fixed=True))
+
+iterN = 10000000
+
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb()]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom()],
+ [Ip2_LudingMat_LudingMat_LudingPhys()],
+ [Law2_ScGeom_LudingPhys_Basic()],
+ ),
+ NewtonIntegrator(damping=0,gravity=[0,0,0]),
+ PyRunner(command='addPlotData()',iterPeriod=1000),
+ PyRunner(command='changeDirection()',iterPeriod=iterN,label="cDir")
+]
+
+velTmp = 0.001
+
+O.bodies[id2].state.vel=[0,0,-velTmp]
+
+def changeDirection():
+ global iterN
+ if (O.bodies[id2].state.vel[2]<0.0):
+ O.bodies[id2].state.vel*=-1.0
+ cDir.iterPeriod = int(iterN/10000.0)
+ elif (O.bodies[id2].state.pos[2]-O.bodies[id1].state.pos[2])-(r1 + r2) > 0.0:
+ iterN = int(iterN*1.2)
+ O.bodies[id2].state.vel[2]*=-1.0
+ cDir.iterPeriod = iterN
+
+def addPlotData():
+ f = [0,0,0]
+ sc = 0
+ try:
+ f=O.forces.f(id2)
+ except:
+ f = [0,0,0]
+
+ s1 = (O.bodies[id2].state.pos[2]-O.bodies[id1].state.pos[2])-(r1 + r2)
+
+ fc1 = f[2]
+ sc1 = -s1/r1
+
+ plot.addData(fc1=fc1, sc=sc1)
+
+plot.plots={'sc':('fc1')}; plot.plot()
+
+from yade import qt
+qt.View()
+
+plot.saveGnuplot('sim-data_LudigPM')
=== added file 'examples/PIDController.py'
--- examples/PIDController.py 1970-01-01 00:00:00 +0000
+++ examples/PIDController.py 2013-07-05 09:11:44 +0000
@@ -0,0 +1,57 @@
+#!/usr/bin/env python
+# encoding: utf-8
+from yade import utils, plot
+
+
+o = Omega()
+fr = 0.0;rho=2000
+tc = 0.001; en = 0.3; et = 0.3; o.dt = 0.02*tc
+
+param = getViscoelasticFromSpheresInteraction(tc,en,et)
+mat1 = O.materials.append(ViscElMat(frictionAngle=fr, density=rho,**param))
+
+spheresID = O.bodies.append(pack.regularHexa(pack.inCylinder((0,0,-2.0),(0,0,2.0),2.0),radius=0.2,gap=0.1,color=(0,1,0),material=mat1))
+
+idWalls = O.bodies.append(geom.facetCylinder(center=(0.0,0.0,0.0),radius = 2.05, height = 4.0, wallMask=6, material=mat1, segmentsNumber = 20, color=(0,0,1)))
+idTop = O.bodies.append(geom.facetCylinder(center=(0.0,0.0,0.0),radius = 2.05, height = 4.0, wallMask=1, material=mat1, segmentsNumber = 5, color=(1,0,0), wire=False))
+
+
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],verletDist=1.0,label='collider'),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
+ [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
+ [Law2_ScGeom_ViscElPhys_Basic()],
+ ),
+ NewtonIntegrator(damping=0,gravity=[0,0,-9.81],label='newtonInt'),
+ TranslationEngine(translationAxis=[0,0,1],velocity=-2.0,ids=idTop,dead=False,label='translat'),
+
+ CombinedKinematicEngine(ids=idTop,label='combEngine',dead=True) +
+ ServoPIDController(axis=[0,0,1],maxVelocity=2.0,iterPeriod=1000,ids=idTop,target=1.0e7,kP=1.0,kI=1.0,kD=1.0) +
+ RotationEngine(rotationAxis=(0,0,1), angularVelocity=10.0, rotateAroundZero=True, zeroPoint=(0,0,0)),
+ PyRunner(command='addPlotData()',iterPeriod=1000, label='graph'),
+ PyRunner(command='switchTranslationEngine()',iterPeriod=45000, nDo = 2, label='switchEng'),
+]
+
+from yade import qt
+qt.View()
+r=qt.Renderer()
+r.bgColor=1,1,1
+
+def addPlotData():
+ fMove = Vector3(0,0,0)
+
+ for i in idTop:
+ fMove += O.forces.f(i)
+
+ plot.addData(z=O.iter, pMove=fMove[2], pFest=fMove[2])
+
+def switchTranslationEngine():
+ print "Switch from TranslationEngine engine to ServoPIDController"
+ translat.dead = True
+ combEngine.dead = False
+
+
+
+plot.plots={'z':('pMove','pFest')}; plot.plot()
=== modified file 'examples/clumps/addToClump-example.py'
--- examples/clumps/addToClump-example.py 2013-03-28 11:01:40 +0000
+++ examples/clumps/addToClump-example.py 2013-07-05 08:21:09 +0000
@@ -48,6 +48,8 @@
keys = b.shape.members.keys()
for ii in range(0,len(keys)):
print '- Body ',keys[ii]
+ print 'inertia:',b.state.inertia
+ print 'mass:',b.state.mass,'\n'
#### show how to use addToClump():
=== modified file 'examples/clumps/releaseFromClump-example.py'
--- examples/clumps/releaseFromClump-example.py 2013-03-28 11:01:40 +0000
+++ examples/clumps/releaseFromClump-example.py 2013-07-11 08:25:03 +0000
@@ -40,9 +40,9 @@
idClump1=O.bodies.clump(bodyList1)
idClump2=O.bodies.clump(bodyList2)
for ii in bodyList1:
- O.bodies[ii].shape.color=(100,100,100)
+ O.bodies[ii].shape.color=(.1,.5,.1)
for ii in bodyList2:
- O.bodies[ii].shape.color=(200,200,200)
+ O.bodies[ii].shape.color=(.2,.2,.7)
#definition for getting informations from all clumps:
def getClumpInfo():
=== modified file 'examples/test/vtk-exporter/vtkExporter.py'
--- examples/test/vtk-exporter/vtkExporter.py 2013-03-28 11:14:53 +0000
+++ examples/test/vtk-exporter/vtkExporter.py 2013-07-10 09:16:15 +0000
@@ -3,12 +3,21 @@
O.bodies.append([
sphere((0,0,0),1),
- sphere((0,2,0),1),
- sphere((0,2,3),2),
+ sphere((0,3,0),1),
+ sphere((0,2,4),2),
+ sphere((0,5,2),1.5),
facet([Vector3(0,-3,-1),Vector3(0,-2,5),Vector3(5,4,0)]),
facet([Vector3(0,-3,-1),Vector3(0,-2,5),Vector3(-5,4,0)])
])
+createInteraction(0,1)
+createInteraction(0,2)
+createInteraction(0,3)
+createInteraction(1,2)
+createInteraction(1,3)
+createInteraction(2,3)
+
vtkExporter = export.VTKExporter('vtkExporterTesting')
vtkExporter.exportSpheres(what=[('dist','b.state.pos.norm()')])
vtkExporter.exportFacets(what=[('pos','b.state.pos')])
+vtkExporter.exportInteractions(what=[('kn','i.phys.kn')])
=== modified file 'examples/triax-tutorial/script-session1.py'
--- examples/triax-tutorial/script-session1.py 2013-04-09 11:46:11 +0000
+++ examples/triax-tutorial/script-session1.py 2013-07-16 09:57:52 +0000
@@ -50,7 +50,7 @@
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
## create walls around the packing
-walls=aabbWalls([mn,mx],thickness=thick,material='walls')
+walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
## use a SpherePack object to generate a random loose particles packing
=== modified file 'lib/base/Math.hpp'
--- lib/base/Math.hpp 2013-06-20 14:56:35 +0000
+++ lib/base/Math.hpp 2013-07-09 12:42:00 +0000
@@ -131,6 +131,19 @@
}
template<typename MatrixT>
+void Matrix_SVD(const MatrixT& in, MatrixT* mU, MatrixT* mS, MatrixT* mV){
+ assert(mU); assert(mS); assert(mV);
+ #if EIGEN_WORLD_VERSION==2 // see Matrix_computeUnitaryPositive
+ Eigen::SVD<MatrixT> svd = Eigen::SVD<MatrixT>(in);
+ #elif EIGEN_WORLD_VERSION==3
+ Eigen::JacobiSVD<MatrixT> svd(in, Eigen::ComputeThinU | Eigen::ComputeThinV);
+ #endif
+ *mU = svd.matrixU();
+ *mV = svd.matrixV();
+ *mS = svd.singularValues().asDiagonal();
+}
+
+template<typename MatrixT>
void matrixEigenDecomposition(const MatrixT& m, MatrixT& mRot, MatrixT& mDiag){
//assert(mRot); assert(mDiag);
Eigen::SelfAdjointEigenSolver<MatrixT> a(m); mRot=a.eigenvectors(); mDiag=a.eigenvalues().asDiagonal();
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2013-06-26 18:15:55 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2013-07-26 18:16:04 +0000
@@ -16,10 +16,6 @@
#include "basicVTKwritter.hpp"
#include "Timer.h"
-#ifdef XVIEW
-#include "Vue3D.h" //FIXME implicit dependencies will look for this class (out of tree) even ifndef XVIEW
-#endif
-
typedef pair<pair<int,int>, vector<double> > Constriction;
using namespace std;
@@ -34,7 +30,7 @@
DECLARE_TESSELATION_TYPES(Network<Tesselation>)
//painfull, but we need that for templates inheritance...
- using _N::T; using _N::x_min; using _N::x_max; using _N::y_min; using _N::y_max; using _N::z_min; using _N::z_max; using _N::Rmoy; using _N::SectionArea; using _N::Height; using _N::Vtotale; using _N::currentTes; using _N::DEBUG_OUT; using _N::nOfSpheres; using _N::x_min_id; using _N::x_max_id; using _N::y_min_id; using _N::y_max_id; using _N::z_min_id; using _N::z_max_id; using _N::boundsIds; using _N::Corner_min; using _N::Corner_max; using _N::Vsolid_tot; using _N::Vtotalissimo; using _N::Vporale; using _N::Ssolid_tot; using _N::V_porale_porosity; using _N::V_totale_porosity; using _N::boundaries; using _N::id_offset; using _N::vtk_infinite_vertices; using _N::vtk_infinite_cells; using _N::num_particles; using _N::fictious_vertex; using _N::boundingCells; using _N::facetVertices;
+ using _N::T; using _N::x_min; using _N::x_max; using _N::y_min; using _N::y_max; using _N::z_min; using _N::z_max; using _N::Rmoy; using _N::SectionArea; using _N::Height; using _N::Vtotale; using _N::currentTes; using _N::DEBUG_OUT; using _N::nOfSpheres; using _N::x_min_id; using _N::x_max_id; using _N::y_min_id; using _N::y_max_id; using _N::z_min_id; using _N::z_max_id; using _N::boundsIds; using _N::Corner_min; using _N::Corner_max; using _N::Vsolid_tot; using _N::Vtotalissimo; using _N::Vporale; using _N::Ssolid_tot; using _N::V_porale_porosity; using _N::V_totale_porosity; using _N::boundaries; using _N::id_offset; using _N::vtk_infinite_vertices; using _N::vtk_infinite_cells; using _N::num_particles; using _N::boundingCells; using _N::facetVertices; using _N::facetNFictious;
//same for functions
using _N::Define_fictious_cells; using _N::AddBoundingPlanes; using _N::boundary;
@@ -42,7 +38,7 @@
FlowBoundingSphere();
bool SLIP_ON_LATERALS;
- bool areaR2Permeability;
+// bool areaR2Permeability;
double TOLERANCE;
double RELAX;
double ks; //Hydraulic Conductivity
@@ -86,7 +82,6 @@
vector <Matrix3r> viscousBodyStress;
vector <Matrix3r> lubBodyStress;
- void mplot ( char *filename);
void Localize();
void Compute_Permeability();
virtual void GaussSeidel (Real dt=0);
@@ -115,7 +110,6 @@
/// Define forces spliting drag and buoyancy terms
void ComputeFacetForcesWithCache(bool onlyCache=false);
void saveVtk ( );
- void MGPost ( );
#ifdef XVIEW
void Dessine_Triangulation ( Vue3D &Vue, RTriangulation &T );
void Dessine_Short_Tesselation ( Vue3D &Vue, Tesselation &Tes );
@@ -151,10 +145,10 @@
void Average_Fluid_Velocity();
void ApplySinusoidalPressure(RTriangulation& Tri, double Amplitude, double Average_Pressure, double load_intervals);
bool isOnSolid (double X, double Y, double Z);
- double MeasurePorePressure (double X, double Y, double Z);
- void MeasurePressureProfile(double Wall_up_y, double Wall_down_y);
- double MeasureAveragedPressure(double Y);
- double MeasureTotalAveragedPressure();
+ double getPorePressure (double X, double Y, double Z);
+ void measurePressureProfile(double Wall_up_y, double Wall_down_y);
+ double averageSlicePressure(double Y);
+ double averagePressure();
double getCell (double X,double Y,double Z);
double boundaryFlux(unsigned int boundaryId);
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2013-07-02 14:45:47 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2013-07-26 18:16:04 +0000
@@ -68,7 +68,6 @@
x_min = 1000.0, x_max = -10000.0, y_min = 1000.0, y_max = -10000.0, z_min = 1000.0, z_max = -10000.0;
currentTes = 0;
nOfSpheres = 0;
- fictious_vertex = 0;
SectionArea = 0, Height=0, Vtotale=0;
vtk_infinite_vertices=0, vtk_infinite_cells=0;
VISCOSITY = 1;
@@ -90,7 +89,7 @@
RAVERAGE = false; /** use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
OUTPUT_BOUDARIES_RADII = false;
RAVERAGE = false; /** if true use the average between the effective radius (inscribed sphere in facet) and the equivalent (circle surface = facet fluid surface) **/
- areaR2Permeability=true;
+// areaR2Permeability=true;
permeability_map = false;
computedOnce=false;
minKdivKmean=0.0001;
@@ -219,9 +218,6 @@
ComputeFacetForcesWithCache();
clock.top("Compute_Forces");
- /** VISUALISATION FILES **/
- MGPost();
-
///*** VUE 3D ***///
#ifdef XVIEW
@@ -406,7 +402,7 @@
return result;
}
template <class Tesselation>
-double FlowBoundingSphere<Tesselation>::MeasurePorePressure (double X, double Y, double Z)
+double FlowBoundingSphere<Tesselation>::getPorePressure (double X, double Y, double Z)
{
if (noCache && T[!currentTes].Max_id()<=0) return 0;//the engine never solved anything
RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
@@ -424,7 +420,7 @@
}
template <class Tesselation>
-void FlowBoundingSphere<Tesselation>::MeasurePressureProfile(double Wall_up_y, double Wall_down_y)
+void FlowBoundingSphere<Tesselation>::measurePressureProfile(double Wall_up_y, double Wall_down_y)
{
if (noCache && T[!currentTes].Max_id()<=0) return;//the engine never solved anything
RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
@@ -450,7 +446,7 @@
}
template <class Tesselation>
-double FlowBoundingSphere<Tesselation>::MeasureAveragedPressure(double Y)
+double FlowBoundingSphere<Tesselation>::averageSlicePressure(double Y)
{
RTriangulation& Tri = T[currentTes].Triangulation();
double P_ave = 0.f;
@@ -468,7 +464,7 @@
return P_ave;
}
template <class Tesselation>
-double FlowBoundingSphere<Tesselation>::MeasureTotalAveragedPressure()
+double FlowBoundingSphere<Tesselation>::averagePressure()
{
RTriangulation& Tri = T[currentTes].Triangulation();
double P = 0.f, Ppond=0.f, Vpond=0.f;
@@ -656,16 +652,25 @@
void FlowBoundingSphere<Tesselation>::Interpolate(Tesselation& Tes, Tesselation& NewTes)
{
Cell_handle old_cell;
-
- RTriangulation& NewTri = NewTes.Triangulation();
RTriangulation& Tri = Tes.Triangulation();
- Finite_cells_iterator cell_end = NewTri.finite_cells_end();
- /*CALCULATION OF VORONOI CENTRES*/
-// if ( !NewTes.Computed() ) NewTes.Compute();
- for (Finite_cells_iterator new_cell = NewTri.finite_cells_begin(); new_cell != cell_end; new_cell++) {
- if (new_cell->info().Pcondition) continue;
- old_cell = Tri.locate((Point) new_cell->info());
- new_cell->info().p() = old_cell->info().p();
+ for (typename Vector_Cell::iterator cell_it=NewTes.cellHandles.begin(); cell_it!=NewTes.cellHandles.end(); cell_it++){
+ Cell_handle& new_cell = *cell_it;
+ if (new_cell->info().Pcondition || new_cell->info().isGhost) continue;
+ Vecteur center ( 0,0,0 );
+ if (new_cell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(new_cell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+ else {
+ Real boundPos=0; int coord=0;
+ for ( int k=0;k<4;k++ ) {
+ if (!new_cell->vertex (k)->info().isFictious) center= center+0.3333333333*(Tes.vertex(new_cell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+ else {
+ coord=boundary (new_cell->vertex(k)->info().id()).coordinate;
+ boundPos=boundary (new_cell->vertex(k)->info().id()).p[coord];
+ }
+ }
+ center=Vecteur(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
+ }
+ old_cell = Tri.locate(Point(center[0],center[1],center[2]));
+ new_cell->info().p() = old_cell->info().shiftedP();
}
// Tes.Clear();//Don't reset to avoid segfault when getting pressure in scripts just after interpolation
}
@@ -703,17 +708,13 @@
int NEG=0, POS=0, pass=0;
bool ref = Tri.finite_cells_begin()->info().isvisited;
- Vecteur n;
-// std::ofstream oFile( "Radii",std::ios::out);
-// std::ofstream fFile( "Radii_Fictious",std::ios::out);
-// std::ofstream kFile ( "LocalPermeabilities" ,std::ios::app );
+// Vecteur n;
Real meanK=0, STDEV=0, meanRadius=0, meanDistance=0;
Real infiniteK=1e10;
- double volume_sub_pore = 0.f;
+// double volume_sub_pore = 0.f;
for (VCell_iterator cell_it=T[currentTes].cellHandles.begin(); cell_it!=T[currentTes].cellHandles.end(); cell_it++){
-// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
Cell_handle& cell = *cell_it;
Point& p1 = cell->info();
for (int j=0; j<4; j++) {
@@ -728,21 +729,16 @@
Sphere& v0 = W[0]->point();
Sphere& v1 = W[1]->point();
Sphere& v2 = W[2]->point();
-#ifdef USE_FAST_MATH
- //FIXME : code not compiling,, do the same as in "else"
- assert((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())>=0 && (W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point())<=1);
- for (int jj=0;jj<3;jj++)
- cell->info().facetSphereCrossSections[j][jj]=0.5*W[jj]->point().weight()*Wm3::FastInvCos1((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point()));
-#else
+
cell->info().facetSphereCrossSections[j]=Vecteur(
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())));
-#endif
+
pass+=1;
Vecteur l = p1 - p2;
distance = sqrt(l.squared_length());
- n = l/distance;
+// n = l/distance;
if (!RAVERAGE) radius = 2* Compute_HydraulicRadius(cell, j);
else radius = (Compute_EffectiveRadius(cell, j)+Compute_EquivalentRadius(cell,j))*0.5;
if (radius<0) NEG++;
@@ -752,7 +748,7 @@
}
// Real h,d;
Real fluidArea=0;
- int test=0;
+// int test=0;
if (distance!=0) {
if (minPermLength>0 && distance_correction) distance=max(minPermLength,distance);
const Vecteur& Surfk = cell->info().facetSurfaces[j];
@@ -777,24 +773,25 @@
meanK += k*k_factor;
if (k<0 && DEBUG_OUT) {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;}
+ 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 "<<endl;}
} else {cout <<"infinite K1!"<<endl; k = infiniteK;}//Will be corrected in the next loop
(cell->info().k_norm())[j]= k*k_factor;
- (neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= k*k_factor;
+// (neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= k*k_factor;
+ if (!neighbour_cell->info().isGhost) (neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]= (cell->info().k_norm())[j];
- if(permeability_map){
- Cell_handle c = cell;
- cell->info().s = cell->info().s + k*distance/fluidArea*this->Volume_Pore_VoronoiFraction (c,j);
- volume_sub_pore += this->Volume_Pore_VoronoiFraction (c,j);}
-
+// if(permeability_map){
+// Cell_handle c = cell;
+// cell->info().s = cell->info().s + k*distance/fluidArea*this->Volume_Pore_VoronoiFraction (c,j);
+// volume_sub_pore += this->Volume_Pore_VoronoiFraction (c,j);}
+//
}
}
cell->info().isvisited = !ref;
- if(permeability_map){cell->info().s = cell->info().s/volume_sub_pore;
- volume_sub_pore = 0.f;}
+// if(permeability_map){cell->info().s = cell->info().s/volume_sub_pore;
+// volume_sub_pore = 0.f;}
}
if (DEBUG_OUT) cout<<"surfneg est "<<surfneg<<endl;
meanK /= pass;
@@ -814,26 +811,21 @@
pass=0;
if (clampKValues) for (VCell_iterator cell_it=T[currentTes].cellHandles.begin(); cell_it!=T[currentTes].cellHandles.end(); cell_it++){
-// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
Cell_handle& cell = *cell_it;
for (int j=0; j<4; j++) {
neighbour_cell = cell->neighbor(j);
if (!Tri.is_infinite(neighbour_cell) && neighbour_cell->info().isvisited==ref) {
pass++;
-// (cell->info().k_norm())[j] = min((cell->info().k_norm())[j], maxKdivKmean*meanK);
(cell->info().k_norm())[j] = max(minKdivKmean*globalK ,min((cell->info().k_norm())[j], maxKdivKmean*globalK));
(neighbour_cell->info().k_norm())[Tri.mirror_index(cell, j)]=(cell->info().k_norm())[j];
-// cout<<(cell->info().k_norm())[j]<<endl;
-// kFile << (cell->info().k_norm())[j] << endl;
}
- }/*cell->info().isvisited = !ref;*/
-
+ }
}
if (DEBUG_OUT) cout << "PassKcorrect = " << pass << endl;
if (DEBUG_OUT) cout << "POS = " << POS << " NEG = " << NEG << " pass = " << pass << endl;
-// A loop to compute the standard deviation of the local K distribution, and use it to include/exclude K values higher then (meanK +/- K_opt_factor*STDEV)
+ // A loop to compute the standard deviation of the local K distribution, and use it to include/exclude K values higher then (meanK +/- K_opt_factor*STDEV)
if (meanKStat)
{
std::ofstream k_opt_file("k_stdev.txt" ,std::ios::out);
@@ -877,21 +869,17 @@
Finite_vertices_iterator vertices_end = Tri.finite_vertices_end();
Real Vgrains = 0;
int grains=0;
-
for (Finite_vertices_iterator V_it = Tri.finite_vertices_begin(); V_it != vertices_end; V_it++) {
- if (!V_it->info().isFictious) {
+ if (!V_it->info().isFictious && !V_it->info().isGhost) {
grains +=1;
- Vgrains += 1.33333333 * M_PI * pow(V_it->point().weight(),1.5);
- }
- }
+ Vgrains += 1.33333333 * M_PI * pow(V_it->point().weight(),1.5);}}
cout<<grains<<"grains - " <<"Vtotale = " << Vtotale << " Vgrains = " << Vgrains << " Vporale1 = " << (Vtotale-Vgrains) << endl;
cout << "Vtotalissimo = " << Vtotalissimo/2 << " Vsolid_tot = " << Vsolid_tot/2 << " Vporale2 = " << Vporale/2 << " Ssolid_tot = " << Ssolid_tot << endl<< endl;
-
if (!RAVERAGE) cout << "------Hydraulic Radius is used for permeability computation------" << endl << endl;
else cout << "------Average Radius is used for permeability computation------" << endl << endl;
cout << "-----Computed_Permeability-----" << endl;}
-// cout << "Negative Permeabilities = " << count_k_neg << endl;
}
+
template <class Tesselation>
vector<double> FlowBoundingSphere<Tesselation>::getConstrictions()
{
@@ -981,14 +969,13 @@
{
RTriangulation& Tri = T[currentTes].Triangulation();
if (Tri.is_infinite(cell->neighbor(j))) return 0;
-
- double Vpore = this->Volume_Pore_VoronoiFraction(cell, j);
- double Ssolid = this->Surface_Solid_Pore(cell, j, SLIP_ON_LATERALS);
+ double Vpore = this->Volume_Pore_VoronoiFraction(cell, j);
+ double Ssolid = this->Surface_Solid_Pore(cell, j, SLIP_ON_LATERALS, /*reuse the same facet data*/ true);
//handle symmetry (tested ok)
- if (SLIP_ON_LATERALS && fictious_vertex>0) {
+ if (SLIP_ON_LATERALS && facetNFictious>0) {
//! Include a multiplier so that permeability will be K/2 or K/4 in symmetry conditions
- Real mult= fictious_vertex==1 ? multSym1 : multSym2;
+ Real mult= facetNFictious==1 ? multSym1 : multSym2;
return Vpore/Ssolid*mult;}
return Vpore/Ssolid;
}
@@ -1087,11 +1074,6 @@
t_dp_max.resize(num_threads);
t_p_max.resize(num_threads);
t_sum_p.resize(num_threads);
-// cerr <<"==== THREADS : "<<num_threads<<" ===="<<endl;
-
-// cout << "tolerance = " << tolerance << endl;
-// cout << "relax = " << relax << endl;
-
Finite_cells_iterator cell_end = Tri.finite_cells_end();
#ifdef GS_OPEN_MP
vector<Finite_cells_iterator> cells_its;
@@ -1380,11 +1362,6 @@
}
vtkfile.end_cells();
-// vtkfile.begin_data("Force",POINT_DATA,VECTORS,FLOAT);
-// for (Finite_vertices_iterator v = T.finite_vertices_begin(); v != T.finite_vertices_end(); ++v)
-// {if (!v->info().isFictious) vtkfile.write_data((v->info().forces)[0],(v->info().forces)[1],(v->info().forces)[2]);}
-// vtkfile.end_data();
-
if (permeability_map){
vtkfile.begin_data("Permeability",CELL_DATA,SCALARS,FLOAT);
for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != Tri.finite_cells_end(); ++cell) {
@@ -1400,7 +1377,6 @@
}
vtkfile.end_data();}
-
if (1){
Average_Relative_Cell_Velocity();
vtkfile.begin_data("Velocity",CELL_DATA,VECTORS,FLOAT);
@@ -1410,44 +1386,7 @@
}
vtkfile.end_data();}
}
-template <class Tesselation>
-void FlowBoundingSphere<Tesselation>::MGPost()
-{
- RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
- Point P;
-
- ofstream file("mgp.out.001");
- file << "<?xml version=\"1.0\"?>" << endl;
- file << "<mgpost mode=\"3D\">" << endl;
- file << " <state time=\"0\">" << endl;
-
- Finite_cells_iterator cell_end = Tri.finite_cells_end();
-
- for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
- if (cell->info().fictious()==0) {
- double p [3] = {0,0,0};
-
- for (int j2=0; j2!=4; j2++) {
- Vertex_handle v = cell->vertex(j2);
- for (int i=0; i<3; i++) {
- p[i] += 0.25* v->point().point()[i];
- }
- }
-
- double pressure = cell->info().p();
- double rad = 0.00025;
-
- file << " <body>" << endl;
- file << " <SPHER r=\"" << rad << "\">" << endl;
- file << " <position x=\"" << p[0] << "\" y=\"" << p[1] << "\" z=\"" << p[2] << "\"/>" << endl;
- file << " <velocity x=\"" << pressure << "\" y=\"" << pressure << "\" z=\"" << pressure << "\"/>" << endl;
- file << " </SPHER>" << endl;
- file << " </body>" << endl;
- }
- }
- file << " </state>" << endl;
- file << "</mgpost>" << endl;
-}
+
#ifdef XVIEW
template <class Tesselation>
void FlowBoundingSphere<Tesselation>::Dessine_Triangulation(Vue3D &Vue, RTriangulation &T)
@@ -1494,26 +1433,8 @@
}
}
}
-template <class Tesselation>
-void FlowBoundingSphere<Tesselation>::mplot (char *filename)
-{
- RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
- std::ofstream plot (filename, std::ios::out);
- Cell_handle permeameter;
- int intervals = 30;
- double Rx = (x_max-x_min) /intervals;
- double Ry = (y_max-y_min) /intervals;
- double Z = (z_max+z_min)/2.0;
- double Press = 0;
- for (double Y=y_min; Y<y_max; Y=Y+Ry) {
- for (double X=x_min; X<x_max; X=X+Rx) {
- permeameter = Tri.locate(Point(X, Y, Z));
- permeameter->info().p()<0? Press=0 : Press=permeameter->info().p();
- plot << Y << " " << X << " " << Press << endl;
- }
- }
-}
-template <class Tesselation>
+
+template <class Tesselation>
double FlowBoundingSphere<Tesselation>::Sample_Permeability(double& x_Min,double& x_Max ,double& y_Min,double& y_Max,double& z_Min,double& z_Max/*, string key*/)
{
double Section = (x_Max-x_Min) * (z_Max-z_Min);
=== modified file 'lib/triangulation/Network.hpp'
--- lib/triangulation/Network.hpp 2013-05-17 17:24:20 +0000
+++ lib/triangulation/Network.hpp 2013-07-26 13:52:11 +0000
@@ -14,17 +14,26 @@
#include "Timer.h"
#include "basicVTKwritter.hpp"
+/**
+Defines class Network. Which contains the geometrical representation of a pore network on the basis of regular triangulation (using CGAL lib)
+The class is the base of the pore-flow model. It has basic functions to compute quantities like void volumes and solid surfaces in the triangulation's elements.
+
+The same data structure is used with different template parameters for periodic and aperiodic boundary conditions. The network is bounded by infinite planes represented in the triangulation by very large spheres (so that their surface looks flat at the scale of the network).
+
+Two triangulations are in fact contained in the network, so that a simulation can switch between them and pass data from one to the other. Otherwise, some info would be lost when the problem is retriangulated.
+*/
+
namespace CGT {
-
+/// Representation of a boundary condition along an axis aligned plane.
struct Boundary
{
- Point p;
- Vecteur normal;
- Vector3r velocity;
- int coordinate;
+ Point p;//position
+ Vecteur normal;//orientation
+ Vector3r velocity;//motion
+ int coordinate;//the axis perpendicular to the boundary
bool flowCondition;//flowCondition=0, pressure is imposed // flowCondition=1, flow is imposed
- Real value;
- bool useMaxMin;
+ Real value;// value of imposed pressure
+ bool useMaxMin;// tells if this boundary was placed following the particles (using min/max of them) or with user defined position
};
@@ -52,17 +61,13 @@
Boundary& boundary (int b) {return boundaries[b-id_offset];}
short id_offset;
int vtk_infinite_vertices, vtk_infinite_cells, num_particles;
-
- int fictious_vertex;
- void AddBoundingPlanes(double altFAR=0);
- void AddBoundingPlane (bool yade, Vecteur Normal, int id_wall, double altFAR=0);
- void AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall, double altFAR=0);
+ void AddBoundingPlanes();
+ void AddBoundingPlane (Vecteur Normal, int id_wall);
+ void AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall );
void Define_fictious_cells( );
- int Detect_facet_fictious_vertices (Cell_handle& cell, int& j);
-
- double Volume_Pore_VoronoiFraction ( Cell_handle& cell, int& j);
+ int detectFacetFictiousVertices (Cell_handle& cell, int& j);
double volumeSolidPore (const Cell_handle& cell);
double volume_single_fictious_pore(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
double volume_double_fictious_pore(const Vertex_handle& SV1, const Vertex_handle& SV2, const Vertex_handle& SV3, const Point& PV1, const Point& PV2, Vecteur& facetSurface);
@@ -72,7 +77,8 @@
Real fast_solid_angle(const Point& STA1, const Point& PTA1, const Point& PTA2, const Point& PTA3);
double volume_double_fictious_pore(Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point PV1);
double volume_single_fictious_pore(Vertex_handle SV1, Vertex_handle SV2, Vertex_handle SV3, Point PV1);
- double Surface_Solid_Pore( Cell_handle cell, int j, bool SLIP_ON_LATERALS);
+ double Volume_Pore_VoronoiFraction ( Cell_handle& cell, int& j, bool reuseFacetData=false);
+ double Surface_Solid_Pore( Cell_handle cell, int j, bool SLIP_ON_LATERALS, bool reuseFacetData=false);
double spherical_triangle_area ( Sphere STA1, Sphere STA2, Sphere STA3, Point PTA1 );
Vecteur surface_double_fictious_facet(Vertex_handle fSV1, Vertex_handle fSV2, Vertex_handle SV3);
@@ -82,9 +88,9 @@
int facetF1, facetF2, facetRe1, facetRe2, facetRe3;
int F1, F2, Re1, Re2;
+ int facetNFictious;
int real_vertex;
- bool facet_detected;
- static const double FAR;
+ double FAR;
static const double ONE_THIRD;
static const int facetVertices [4][3];
static const int permut3 [3][3];
=== modified file 'lib/triangulation/Network.ipp'
--- lib/triangulation/Network.ipp 2013-05-17 17:24:20 +0000
+++ lib/triangulation/Network.ipp 2013-07-26 13:52:11 +0000
@@ -19,13 +19,7 @@
// using namespace boost;
namespace CGT {
-// template<class Tesselation> const double Network<Tesselation>::FAR = 50000;
-// template<class Tesselation> const double Network<Tesselation>::ONE_THIRD = 1.0/3.0;
-// template<class Tesselation> const int Network<Tesselation>::facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
-// template<class Tesselation> const int Network<Tesselation>::permut3 [3][3] = {{0,1,2},{1,2,0},{2,0,1}};
-// template<class Tesselation> const int Network<Tesselation>::permut4 [4][4] = {{0,1,2,3},{1,2,3,0},{2,3,0,1},{3,0,1,2}};
-
- template<class Tesselation> const double Network<Tesselation>::FAR = 50000;
+// template<class Tesselation> const double Network<Tesselation>::FAR = 50000;
template<class Tesselation> const double Network<Tesselation>::ONE_THIRD = 1.0/3.0;
template<class Tesselation> const int Network<Tesselation>::facetVertices [4][3] = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
template<class Tesselation> const int Network<Tesselation>::permut3 [3][3] = {{0,1,2},{1,2,0},{2,0,1}};
@@ -36,16 +30,15 @@
template<class Tesselation>
Network<Tesselation>::Network(){
+ FAR = 50000;
facetF1=facetF2=facetRe1=facetRe2=facetRe3=0;
F1=F2=Re1=Re2=0;
- facet_detected = false;
-
}
template<class Tesselation>
-int Network<Tesselation>::Detect_facet_fictious_vertices (Cell_handle& cell, int& j)
+int Network<Tesselation>::detectFacetFictiousVertices (Cell_handle& cell, int& j)
{
- fictious_vertex = 0;
+ facetNFictious = 0;
int real_vertex=0;
Sphere v [3];
Vertex_handle W [3];
@@ -54,28 +47,27 @@
W[kk] = cell->vertex(facetVertices[j][kk]);
v[kk] = cell->vertex(facetVertices[j][kk])->point();
if (W[kk]->info().isFictious) {
- if (fictious_vertex==0) {F1=facetVertices[j][kk];facetF1=kk;} else
+ if (facetNFictious==0) {F1=facetVertices[j][kk];facetF1=kk;} else
{F2 = facetVertices[j][kk];facetF2=kk;}
- fictious_vertex +=1;
+ facetNFictious +=1;
} else {
if (real_vertex==0) {Re1=facetVertices[j][kk];facetRe1=kk;} else if (real_vertex==1)
{Re2=facetVertices[j][kk];facetRe2=kk;} else if (real_vertex==2)
{Re2=facetVertices[j][kk];facetRe3=kk;}
real_vertex+=1;}}
- facet_detected = true;
- return fictious_vertex;
+ return facetNFictious;
}
template<class Tesselation>
-double Network<Tesselation>::Volume_Pore_VoronoiFraction (Cell_handle& cell, int& j)
+double Network<Tesselation>::Volume_Pore_VoronoiFraction (Cell_handle& cell, int& j, bool reuseFacetData)
{
Point& p1 = cell->info();
Point& p2 = cell->neighbor(j)->info();
- fictious_vertex = Detect_facet_fictious_vertices (cell,j);
+ if (!reuseFacetData) facetNFictious = detectFacetFictiousVertices (cell,j);
Sphere v [3];
Vertex_handle W [3];
for (int kk=0; kk<3; kk++) {W[kk] = cell->vertex(facetVertices[j][kk]);v[kk] = cell->vertex(facetVertices[j][kk])->point();}
- switch (fictious_vertex) {
+ switch (facetNFictious) {
case (0) : {
Vertex_handle& SV1 = W[0];
Vertex_handle& SV2 = W[1];
@@ -269,11 +261,9 @@
}
template<class Tesselation>
-double Network<Tesselation>::Surface_Solid_Pore(Cell_handle cell, int j, bool SLIP_ON_LATERALS)
+double Network<Tesselation>::Surface_Solid_Pore(Cell_handle cell, int j, bool SLIP_ON_LATERALS, bool reuseFacetData)
{
- if (!facet_detected) fictious_vertex=Detect_facet_fictious_vertices(cell,j);
-
-// RTriangulation& Tri = T[currentTes].Triangulation();
+ if (!reuseFacetData) facetNFictious=detectFacetFictiousVertices(cell,j);
Point& p1 = cell->info();
Point& p2 = cell->neighbor(j)->info();
@@ -287,7 +277,7 @@
W[kk] = cell->vertex(facetVertices[j][kk]);
v[kk] = cell->vertex(facetVertices[j][kk])->point();}
- switch (fictious_vertex) {
+ switch (facetNFictious) {
case (0) : {
Vertex_handle& SV1 = W[0];
Vertex_handle& SV2 = W[1];
@@ -449,10 +439,9 @@
}
template<class Tesselation>
-void Network<Tesselation>::AddBoundingPlanes(double altFAR)
+void Network<Tesselation>::AddBoundingPlanes()
{
Tesselation& Tes = T[currentTes];
-
//FIXME: Id's order in boundsIds is done according to the enumerotation of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
y_min_id = Tes.Max_id() + 2;
boundsIds[0]=&y_min_id;
@@ -472,54 +461,44 @@
id_offset = Tes.Max_id() +1;//so that boundaries[vertex->id - offset] gives the ordered boundaries (also see function Boundary& boundary(int b))
- AddBoundingPlane (true, Vecteur(0,1,0) , y_min_id, altFAR);
- AddBoundingPlane (true, Vecteur(0,-1,0) , y_max_id, altFAR);
- AddBoundingPlane (true, Vecteur(-1,0,0) , x_max_id, altFAR);
- AddBoundingPlane (true, Vecteur(1,0,0) , x_min_id, altFAR);
- AddBoundingPlane (true, Vecteur(0,0,1) , z_min_id, altFAR);
- AddBoundingPlane (true, Vecteur(0,0,-1) , z_max_id, altFAR);
+ AddBoundingPlane (Vecteur(0,1,0) , y_min_id);
+ AddBoundingPlane (Vecteur(0,-1,0) , y_max_id);
+ AddBoundingPlane (Vecteur(-1,0,0) , x_max_id);
+ AddBoundingPlane (Vecteur(1,0,0) , x_min_id);
+ AddBoundingPlane (Vecteur(0,0,1) , z_min_id);
+ AddBoundingPlane (Vecteur(0,0,-1) , z_max_id);
// AddBoundingPlanes(true);
}
template<class Tesselation>
-void Network<Tesselation>::AddBoundingPlane (bool yade, Vecteur Normal, int id_wall, double altFAR)
+void Network<Tesselation>::AddBoundingPlane (Vecteur Normal, int id_wall)
{
- Tesselation& Tes = T[currentTes];
- if (altFAR==0) altFAR=FAR;
+// Tesselation& Tes = T[currentTes];
+ //FIXME: pre-condition: the normal is axis-aligned
int Coordinate = abs(Normal[0])*0 + abs(Normal[1])*1 + abs(Normal[2])*2;
double pivot = Normal[Coordinate]<0 ?
Corner_max.x()*abs(Normal[0])+Corner_max.y()*abs(Normal[1])+Corner_max.z()*abs(Normal[2]) : Corner_min.x()*abs(Normal[0])+Corner_min.y()*abs(Normal[1])+Corner_min.z()*abs(Normal[2]);
-
- Tes.insert(0.5*(Corner_min.x() +Corner_max.x())*(1-abs(Normal[0]))+(pivot-Normal[0]*altFAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]),
- 0.5*(Corner_max.y() +Corner_min.y())*(1-abs(Normal[1]))+(pivot-Normal[1]*altFAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]),
- 0.5*(Corner_max.z() +Corner_min.z())*(1-abs(Normal[2]))+(pivot-Normal[2]*altFAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]),
- altFAR*(Corner_max.y()-Corner_min.y()), id_wall, true);
- if (Normal[Coordinate]<0) boundaries[id_wall-id_offset].p = Corner_max;
- else boundaries[id_wall-id_offset].p = Corner_min;
-
- boundaries[id_wall-id_offset].normal = Normal;
- boundaries[id_wall-id_offset].coordinate = Coordinate;
-
- boundaries[id_wall-id_offset].flowCondition = 1;
- boundaries[id_wall-id_offset].value = 0;
-
- if(DEBUG_OUT) cout << "A boundary -max/min-has been created. ID = " << id_wall << " position = " << 0.5*(Corner_min.x() +Corner_max.x())*(1-abs(Normal[0]))+(pivot-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]) << " , " << 0.5*(Corner_max.y() +Corner_min.y())*(1-abs(Normal[1]))+(pivot-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]) << " , " << 0.5*(Corner_max.z() +Corner_min.z())*(1-abs(Normal[2]))+(pivot-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]) << ". Radius = " << FAR*(Corner_max.y()-Corner_min.y()) << endl;
+ Real center [3] ={ 0.5*(Corner_min.x() +Corner_max.x())*(1-abs(Normal[0]))+pivot*abs(Normal[0]),
+ 0.5*(Corner_max.y() +Corner_min.y())*(1-abs(Normal[1]))+pivot*abs(Normal[1]),
+ 0.5*(Corner_max.z() +Corner_min.z())*(1-abs(Normal[2]))+pivot*abs(Normal[2])};
+
+ AddBoundingPlane(center,0,Normal,id_wall);
}
template<class Tesselation>
-void Network<Tesselation>::AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall, double altFAR )
+void Network<Tesselation>::AddBoundingPlane (Real center[3], double thickness, Vecteur Normal, int id_wall )
{
Tesselation& Tes = T[currentTes];
- if (altFAR==0) altFAR=FAR;
+
int Coordinate = abs(Normal[0])*0 + abs(Normal[1])*1 + abs(Normal[2])*2;
- Tes.insert((center[0]+Normal[0]*thickness/2)*(1-abs(Normal[0])) + (center[0]+Normal[0]*thickness/2-Normal[0]*altFAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]),
- (center[1]+Normal[1]*thickness/2)*(1-abs(Normal[1])) + (center[1]+Normal[1]*thickness/2-Normal[1]*altFAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]),
- (center[2]+Normal[2]*thickness/2)*(1-abs(Normal[2])) + (center[2]+Normal[2]*thickness/2-Normal[2]*altFAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]),
- altFAR*(Corner_max.y()-Corner_min.y()), id_wall, true);
+ Tes.insert((center[0]+Normal[0]*thickness/2)*(1-abs(Normal[0])) + (center[0]+Normal[0]*thickness/2-Normal[0]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[0]),
+ (center[1]+Normal[1]*thickness/2)*(1-abs(Normal[1])) + (center[1]+Normal[1]*thickness/2-Normal[1]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[1]),
+ (center[2]+Normal[2]*thickness/2)*(1-abs(Normal[2])) + (center[2]+Normal[2]*thickness/2-Normal[2]*FAR*(Corner_max.y()-Corner_min.y()))*abs(Normal[2]),
+ FAR*(Corner_max.y()-Corner_min.y()), id_wall, true);
Point P (center[0],center[1],center[2]);
boundaries[id_wall-id_offset].p = P;
=== modified file 'lib/triangulation/def_types.h'
--- lib/triangulation/def_types.h 2013-07-26 14:55:13 +0000
+++ lib/triangulation/def_types.h 2013-07-30 17:03:30 +0000
@@ -8,7 +8,6 @@
#include <yade/lib/base/Math.hpp>
//Define basic types from CGAL templates
-
#ifndef _Def_types
#define _Def_types
@@ -23,7 +22,6 @@
#include <CGAL/number_utils.h>
#include <boost/static_assert.hpp>
-
namespace CGT {
//Robust kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
@@ -132,6 +130,7 @@
pression=0;
invVoidV=0;
fict=0;
+ isGhost=false;
}
bool isGhost;
double inv_sum_k;
@@ -171,6 +170,8 @@
FlowVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
FlowVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
Vecteur forces;
+ bool isGhost;
+ FlowVertexInfo (void) {isGhost=false;}
inline Vecteur& force (void) {return forces;}
inline Vecteur& vel (void) {return Grain_Velocity;}
inline Real& vol_cells (void) {return volume_incident_cells;}
@@ -192,7 +193,6 @@
PeriodicCellInfo (void){
_pression=&pression;
period[0]=period[1]=period[2]=0;
- isGhost=false;
baseIndex=-1;
volumeSign=0;}
~PeriodicCellInfo (void) {}
@@ -211,7 +211,6 @@
PeriodicVertexInfo& operator= (const Vecteur &u) { Vecteur::operator= (u); return *this; }
PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
PeriodicVertexInfo& operator= (const unsigned int &id) { i= id; return *this; }
- bool isGhost;
int period[3];
//FIXME: the name is misleading, even non-ghost can be out of the period and therefore they need to be shifted as well
inline const Vecteur ghostShift (void) {
=== modified file 'pkg/common/KinematicEngines.cpp'
--- pkg/common/KinematicEngines.cpp 2013-03-07 18:28:01 +0000
+++ pkg/common/KinematicEngines.cpp 2013-07-04 17:36:34 +0000
@@ -4,7 +4,7 @@
#include<yade/pkg/dem/Shop.hpp>
#include<yade/lib/smoothing/LinearInterpolate.hpp>
-YADE_PLUGIN((KinematicEngine)(CombinedKinematicEngine)(TranslationEngine)(HarmonicMotionEngine)(RotationEngine)(HelixEngine)(InterpolatingHelixEngine)(HarmonicRotationEngine));
+YADE_PLUGIN((KinematicEngine)(CombinedKinematicEngine)(TranslationEngine)(HarmonicMotionEngine)(RotationEngine)(HelixEngine)(InterpolatingHelixEngine)(HarmonicRotationEngine)(ServoPIDController));
CREATE_LOGGER(KinematicEngine);
@@ -142,3 +142,43 @@
RotationEngine::apply(ids);
}
+void ServoPIDController::apply(const vector<Body::id_t>& ids){
+
+ if (iterPrevStart < 0 or ((scene->iter-iterPrevStart)>=iterPeriod)) {
+
+ Vector3r tmpForce = Vector3r::Zero();
+
+ if (ids.size()>0) {
+ FOREACH(Body::id_t id,ids){
+ assert(id<(Body::id_t)scene->bodies->size());
+ tmpForce += scene->forces.getForce(id);
+ }
+ } else {
+ LOG_WARN("The list of ids is empty!");
+ }
+
+ axis.normalize();
+ tmpForce = tmpForce.cwiseProduct(axis); // Take into account given axis
+ errorCur = tmpForce.norm() - target; // Find error
+
+ const Real pTerm = errorCur*kP; // Calculate proportional term
+ iTerm += errorCur*kI; // Calculate integral term
+ const Real dTerm = (errorCur-errorPrev)*kD; // Calculate derivative term
+
+ errorPrev = errorCur; // Save the current value of the error
+
+ curVel = (pTerm + iTerm + dTerm); // Calculate current velocity
+
+ if (fabs(curVel) > fabs(maxVelocity)) {
+ curVel*=fabs(maxVelocity)/fabs(curVel);
+ }
+
+ iterPrevStart = scene->iter;
+ }
+
+ translationAxis = axis;
+ velocity = curVel;
+
+ TranslationEngine::apply(ids);
+}
+
=== modified file 'pkg/common/KinematicEngines.hpp'
--- pkg/common/KinematicEngines.hpp 2012-09-21 19:03:18 +0000
+++ pkg/common/KinematicEngines.hpp 2013-07-04 17:36:34 +0000
@@ -93,3 +93,25 @@
};
REGISTER_SERIALIZABLE(HarmonicRotationEngine);
+struct ServoPIDController: public TranslationEngine{
+ virtual void apply(const vector<Body::id_t>& ids);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(ServoPIDController,TranslationEngine,"PIDController servo-engine for applying prescribed force on bodies. http://en.wikipedia.org/wiki/PID_controller";,
+ ((Real,maxVelocity,0.0,,"Velocity [m/s]"))
+ ((Vector3r,axis,Vector3r::Zero(),,"Unit vector along which apply the velocity [-]"))
+ ((Real,target,0.0,,"Target value for the controller [N]"))
+ ((Real,kP,0.0,,"Proportional gain/coefficient for the PID-controller [-]"))
+ ((Real,kI,0.0,,"Integral gain/coefficient for the PID-controller [-]"))
+ ((Real,kD,0.0,,"Derivative gain/coefficient for the PID-controller [-]"))
+ ((Real,iTerm,0.0,,"Integral term [N]"))
+ ((Real,curVel,0.0,,"Current applied velocity [m/s]"))
+ ((Real,errorCur,0.0,,"Current error [N]"))
+ ((Real,errorPrev,0.0,,"Previous error [N]"))
+ ((long,iterPeriod,100.0,,"Periodicity criterion of velocity correlation [-]"))
+ ((long,iterPrevStart,-1.0,,"Previous iteration of velocity correlation [-]"))
+ /* attrs*/,
+ /* ctor */,
+ /* py */
+ )
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(ServoPIDController);
=== modified file 'pkg/dem/BubbleMat.hpp'
--- pkg/dem/BubbleMat.hpp 2013-06-15 19:17:48 +0000
+++ pkg/dem/BubbleMat.hpp 2013-07-05 15:54:47 +0000
@@ -11,8 +11,8 @@
/********************** BubbleMat ****************************/
class BubbleMat : public Material {
public:
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(BubbleMat,Material,"TODO DOC",
- ((Real,surfaceTension,1/*TODO some realistic value*/,,"TODO DOC"))
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(BubbleMat,Material,"material for bubble interactions, for use with other Bubble classes",
+ ((Real,surfaceTension,0.07197,,"The surface tension in the fluid surrounding the bubbles. The default value is that of water at 25 degrees Celcius."))
,
createIndex();
density=1; // TODO density default value
@@ -29,17 +29,17 @@
static Real computeForce(Real penetrationDepth, Real surfaceTension, Real rAvg, int newtonIter, Real newtonTol);
virtual ~BubblePhys(){};
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(BubblePhys,IPhys,"TODO DOC",
- ((Vector3r,normalForce,Vector3r::Zero(),,"TODO DOC"))
- ((Real,surfaceTension,NaN,,"TODO DOC"))
- ((Real,fN,NaN,,"TODO DOC"))
- ((Real,rAvg,NaN,,"TODO DOC"))
- ((Real,newtonIter,50,,"TODO DOC"))
- ((Real,newtonTol,1e-6,,"TODO DOC"))
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(BubblePhys,IPhys,"Physics of bubble-bubble interactions, for use with BubbleMat",
+ ((Vector3r,normalForce,Vector3r::Zero(),,"Normal force"))
+ ((Real,surfaceTension,NaN,,"Surface tension of the surrounding liquid"))
+ ((Real,fN,NaN,,"Contact normal force"))
+ ((Real,rAvg,NaN,,"Average radius of the two interacting bubbles"))
+ ((Real,newtonIter,50,,"Maximum number of force iterations allowed"))
+ ((Real,newtonTol,1e-6,,"Convergence criteria for force iterations"))
,
createIndex();
,
- .def("computeForce",&BubblePhys::computeForce,"TODO DOC")
+ .def("computeForce",&BubblePhys::computeForce,"Computes the normal force acting between the two interacting bubbles using the Newton-Rhapson method")
.staticmethod("computeForce")
);
DECLARE_LOGGER;
@@ -55,7 +55,7 @@
virtual void go(const shared_ptr<Material>& m1, const shared_ptr<Material>& m2, const shared_ptr<Interaction>& interaction);
FUNCTOR2D(BubbleMat,BubbleMat);
DECLARE_LOGGER;
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_BubbleMat_BubbleMat_BubblePhys,IPhysFunctor,"TODO DOC",
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_BubbleMat_BubbleMat_BubblePhys,IPhysFunctor,"Generates bubble interactions.Used in the contact law Law2_ScGeom_BubblePhys_Bubble.",
);
};
REGISTER_SERIALIZABLE(Ip2_BubbleMat_BubbleMat_BubblePhys);
@@ -66,7 +66,7 @@
public:
void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* interaction);
FUNCTOR2D(GenericSpheresContact,BubblePhys);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_BubblePhys_Bubble,LawFunctor,"TODO DOC",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_BubblePhys_Bubble,LawFunctor,"Constitutive law for Bubble model.",
,
/*ctor*/,
);
=== removed file 'pkg/dem/CohesiveStateRPMRecorder.cpp'
--- pkg/dem/CohesiveStateRPMRecorder.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/CohesiveStateRPMRecorder.cpp 1970-01-01 00:00:00 +0000
@@ -1,19 +0,0 @@
-#include"CohesiveStateRPMRecorder.hpp"
-
-YADE_PLUGIN((CohesiveStateRPMRecorder));
-CREATE_LOGGER(CohesiveStateRPMRecorder);
-
-void CohesiveStateRPMRecorder::action() {
- numberCohesiveContacts=0;
- //Check all interactions
- FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){
- if(!i->isReal()) continue; //Check whether they are real
- const shared_ptr<RpmPhys>& contPhys = YADE_PTR_CAST<RpmPhys>(i->phys);
- if (contPhys->isCohesive==true) { //Check whether they are cohesive
- numberCohesiveContacts++; //If yes - calculate them
- }
- }
- //Save data to a file
- out<<scene->iter<<" "<<numberCohesiveContacts<<"\n";
- out.close();
-}
=== removed file 'pkg/dem/CohesiveStateRPMRecorder.hpp'
--- pkg/dem/CohesiveStateRPMRecorder.hpp 2010-11-12 08:03:16 +0000
+++ pkg/dem/CohesiveStateRPMRecorder.hpp 1970-01-01 00:00:00 +0000
@@ -1,19 +0,0 @@
-/*! This class is for recording number of cohesive contacts
- * of RPM model in a file. Class derived from Recorder
-*/
-#pragma once
-#include <yade/pkg/common/Recorder.hpp>
-#include <yade/pkg/common/Dispatching.hpp>
-#include <yade/pkg/dem/RockPM.hpp>
-
-class CohesiveStateRPMRecorder: public Recorder {
- std::ofstream outFile;
- public:
- virtual void action();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohesiveStateRPMRecorder,Recorder,
- "Store number of cohesive contacts in RPM model to file.",
- ((int,numberCohesiveContacts,0,,"Number of cohesive contacts found at last run. [-]")),
- initRun=true;);
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(CohesiveStateRPMRecorder);
=== modified file 'pkg/dem/FlowEngine.cpp'
--- pkg/dem/FlowEngine.cpp 2013-05-06 15:50:31 +0000
+++ pkg/dem/FlowEngine.cpp 2013-07-26 18:16:04 +0000
@@ -58,13 +58,13 @@
timingDeltas->checkpoint ( "Update_Volumes" );
Eps_Vol_Cumulative += eps_vol_max;
- if ( (EpsVolPercent_RTRG>0 && Eps_Vol_Cumulative > EpsVolPercent_RTRG) || retriangulationLastIter>PermuteInterval) {
- Update_Triangulation = true;
+ if ( (defTolerance>0 && Eps_Vol_Cumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval) {
+ updateTriangulation = true;
Eps_Vol_Cumulative=0;
retriangulationLastIter=0;
ReTrg++;
} else {
- Update_Triangulation = false;
+ updateTriangulation = false;
retriangulationLastIter++;}
@@ -95,11 +95,11 @@
timingDeltas->checkpoint ( "Applying Forces" );
int sleeping = 0;
if (multithread && !first) {
- while (Update_Triangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/
+ while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/
sleeping++;
boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
if (Debug && sleeping) cerr<<"sleeping..."<<sleeping<<endl;
- if (Update_Triangulation || (ellapsedIter>(0.5*PermuteInterval) && backgroundCompleted)) {
+ if (updateTriangulation || (ellapsedIter>(0.5*meshUpdateInterval) && backgroundCompleted)) {
if (Debug) cerr<<"switch flow solver"<<endl;
if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
if (fluidBulkModulus>0) solver->Interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
@@ -112,7 +112,7 @@
setPositionsBuffer(false);//set "parallel" buffer for background calculation
backgroundCompleted=false;
retriangulationLastIter=ellapsedIter;
- Update_Triangulation=false;
+ updateTriangulation=false;
ellapsedIter=0;
boost::thread workerThread(&FlowEngine::backgroundAction,this);
workerThread.detach();
@@ -126,13 +126,12 @@
ellapsedIter++;
}
} else {
- if (Update_Triangulation && !first) {
+ if (updateTriangulation && !first) {
Build_Triangulation (P_zero, solver);
Initialize_volumes(solver);
ComputeViscousForces(*solver);
- Update_Triangulation = false;}
+ updateTriangulation = false;}
}
- if ( velocity_profile ) /*flow->FluidVelocityProfile();*/solver->Average_Fluid_Velocity();
first=false;
timingDeltas->checkpoint ( "Triangulate + init volumes" );
}
@@ -153,36 +152,12 @@
void FlowEngine::BoundaryConditions ( Solver& flow )
{
- if ( flow->y_min_id>=0 ) {
- flow->boundary ( flow->y_min_id ).flowCondition=Flow_imposed_BOTTOM_Boundary;
- flow->boundary ( flow->y_min_id ).value=Pressure_BOTTOM_Boundary;
- flow->boundary ( flow->y_min_id ).velocity = Vector3r::Zero();
- }
- if ( flow->y_max_id>=0 ) {
- flow->boundary ( flow->y_max_id ).flowCondition=Flow_imposed_TOP_Boundary;
- flow->boundary ( flow->y_max_id ).value=Pressure_TOP_Boundary;
- flow->boundary ( flow->y_max_id ).velocity = topBoundaryVelocity;
- }
- if ( flow->x_max_id>=0 ) {
- flow->boundary ( flow->x_max_id ).flowCondition=Flow_imposed_RIGHT_Boundary;
- flow->boundary ( flow->x_max_id ).value=Pressure_RIGHT_Boundary;
- flow->boundary ( flow->x_max_id ).velocity = Vector3r::Zero();
- }
- if ( flow->x_min_id>=0 ) {
- flow->boundary ( flow->x_min_id ).flowCondition=Flow_imposed_LEFT_Boundary;
- flow->boundary ( flow->x_min_id ).value=Pressure_LEFT_Boundary;
- flow->boundary ( flow->x_min_id ).velocity = Vector3r::Zero();
- }
- if ( flow->z_max_id>=0 ) {
- flow->boundary ( flow->z_max_id ).flowCondition=Flow_imposed_FRONT_Boundary;
- flow->boundary ( flow->z_max_id ).value=Pressure_FRONT_Boundary;
- flow->boundary ( flow->z_max_id ).velocity = Vector3r::Zero();
- }
- if ( flow->z_min_id>=0 ) {
- flow->boundary ( flow->z_min_id ).flowCondition=Flow_imposed_BACK_Boundary;
- flow->boundary ( flow->z_min_id ).value=Pressure_BACK_Boundary;
- flow->boundary ( flow->z_min_id ).velocity = Vector3r::Zero();
- }
+
+ for (int k=0;k<6;k++) {
+ flow->boundary (wallIds[k]).flowCondition=!bndCondIsPressure[k];
+ flow->boundary (wallIds[k]).value=bndCondValue[k];
+ flow->boundary (wallIds[k]).velocity = boundaryVelocity[k];//FIXME: needs correct implementation, maybe update the cached pos/vel?
+ }
}
template<class Solver>
@@ -222,7 +197,7 @@
{
flow->Vtotalissimo=0; flow->Vsolid_tot=0; flow->Vporale=0; flow->Ssolid_tot=0;
flow->SLIP_ON_LATERALS=slip_boundary;
- flow->k_factor = permeability_factor;
+ flow->k_factor = permeabilityFactor;
flow->DEBUG_OUT = Debug;
flow->useSolver = useSolver;
#ifdef EIGENSPARSE_LIB
@@ -231,7 +206,6 @@
#endif
flow->meanKStat = meanKStat;
flow->VISCOSITY = viscosity;
- flow->areaR2Permeability=areaR2Permeability;
flow->TOLERANCE=Tolerance;
flow->RELAX=Relax;
flow->clampKValues = clampKValues;
@@ -276,18 +250,14 @@
for ( Finite_cells_iterator cell = flow->T[flow->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ )
flow->T[flow->currentTes].cellHandles.push_back(cell);
flow->DisplayStatistics ();
- flow->Compute_Permeability ( );
-
+ flow->Compute_Permeability();
porosity = flow->V_porale_porosity/flow->V_totale_porosity;
- if ( compute_K ) {BoundaryConditions ( flow ); K = flow->Sample_Permeability ( flow->x_min, flow->x_max, flow->y_min, flow->y_max, flow->z_min, flow->z_max );}
-
BoundaryConditions ( flow );
flow->Initialize_pressures ( P_zero );
if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
- if ( WaveAction ) flow->ApplySinusoidalPressure ( flow->T[flow->currentTes].Triangulation(), Sinus_Amplitude, Sinus_Average, 30 );
-
+ if ( WaveAction ) flow->ApplySinusoidalPressure ( flow->T[flow->currentTes].Triangulation(), sineMagnitude, sineAverage, 30 );
if ( viscousShear || normalLubrication || shearLubrication) flow->computeEdgesSurfaces();
}
@@ -334,19 +304,12 @@
flow->id_offset = id_offset;
flow->SectionArea = ( flow->x_max - flow->x_min ) * ( flow->z_max-flow->z_min );
flow->Vtotale = ( flow->x_max-flow->x_min ) * ( flow->y_max-flow->y_min ) * ( flow->z_max-flow->z_min );
- flow->y_min_id=wallBottomId;
- flow->y_max_id=wallTopId;
- flow->x_max_id=wallRightId;
- flow->x_min_id=wallLeftId;
- flow->z_min_id=wallBackId;
- flow->z_max_id=wallFrontId;
-
- if ( flow->y_min_id>=0 ) flow->boundary ( flow->y_min_id ).useMaxMin = BOTTOM_Boundary_MaxMin;
- if ( flow->y_max_id>=0 ) flow->boundary ( flow->y_max_id ).useMaxMin = TOP_Boundary_MaxMin;
- if ( flow->x_max_id>=0 ) flow->boundary ( flow->x_max_id ).useMaxMin = RIGHT_Boundary_MaxMin;
- if ( flow->x_min_id>=0 ) flow->boundary ( flow->x_min_id ).useMaxMin = LEFT_Boundary_MaxMin;
- if ( flow->z_max_id>=0 ) flow->boundary ( flow->z_max_id ).useMaxMin = FRONT_Boundary_MaxMin;
- if ( flow->z_min_id>=0 ) flow->boundary ( flow->z_min_id ).useMaxMin = BACK_Boundary_MaxMin;
+ flow->y_min_id=wallIds[ymin];
+ flow->y_max_id=wallIds[ymax];
+ flow->x_max_id=wallIds[xmax];
+ flow->x_min_id=wallIds[xmin];
+ flow->z_min_id=wallIds[zmin];
+ flow->z_max_id=wallIds[zmax];
//FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
flow->boundsIds[0]= &flow->x_min_id;
@@ -356,20 +319,18 @@
flow->boundsIds[4]= &flow->z_min_id;
flow->boundsIds[5]= &flow->z_max_id;
+ for (int k=0;k<6;k++) flow->boundary ( *flow->boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k];
+
+// if ( flow->y_min_id>=0 ) flow->boundary ( flow->y_min_id ).useMaxMin = boundaryUseMaxMin[ymin];
+// if ( flow->y_max_id>=0 ) flow->boundary ( flow->y_max_id ).useMaxMin = boundaryUseMaxMin[ymax];
+// if ( flow->x_max_id>=0 ) flow->boundary ( flow->x_max_id ).useMaxMin = boundaryUseMaxMin[xmax];
+// if ( flow->x_min_id>=0 ) flow->boundary ( flow->x_min_id ).useMaxMin = boundaryUseMaxMin[xmin];
+// if ( flow->z_max_id>=0 ) flow->boundary ( flow->z_max_id ).useMaxMin = boundaryUseMaxMin[zmax];
+// if ( flow->z_min_id>=0 ) flow->boundary ( flow->z_min_id ).useMaxMin = boundaryUseMaxMin[zmin];
+
flow->Corner_min = CGT::Point ( flow->x_min, flow->y_min, flow->z_min );
flow->Corner_max = CGT::Point ( flow->x_max, flow->y_max, flow->z_max );
-
- if ( Debug ) {
- cout << "Section area = " << flow->SectionArea << endl;
- cout << "Vtotale = " << flow->Vtotale << endl;
- cout << "x_min = " << flow->x_min << endl;
- cout << "x_max = " << flow->x_max << endl;
- cout << "y_max = " << flow->y_max << endl;
- cout << "y_min = " << flow->y_min << endl;
- cout << "z_min = " << flow->z_min << endl;
- cout << "z_max = " << flow->z_max << endl;
- cout << endl << "Adding Boundary------" << endl;
- }
+
//assign BCs types and values
BoundaryConditions ( flow );
@@ -377,9 +338,10 @@
for ( int i=0; i<6; i++ ) {
if ( *flow->boundsIds[i]<0 ) continue;
CGT::Vecteur Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
- if ( flow->boundary ( *flow->boundsIds[i] ).useMaxMin ) flow->AddBoundingPlane ( true, Normal, *flow->boundsIds[i] );
+ if ( flow->boundary ( *flow->boundsIds[i] ).useMaxMin ) flow->AddBoundingPlane(Normal, *flow->boundsIds[i] );
else {
for ( int h=0;h<3;h++ ) center[h] = buffer[*flow->boundsIds[i]].pos[h];
+// cerr << "id="<<*flow->boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl;
flow->AddBoundingPlane ( center, wall_thickness, Normal,*flow->boundsIds[i] );
}
}
@@ -431,7 +393,6 @@
case ( 3 ) : cell->info().volume() = Volume_cell_triple_fictious ( cell ); break;
default: break;
}
-
if (flow->fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow->volumeSolidPore(cell) ); }
}
if (Debug) cout << "Volumes initialised." << endl;
@@ -493,13 +454,13 @@
dVol=cell->info().volumeSign* ( newVol - cell->info().volume() );
cell->info().dv() = dVol*invDeltaT;
cell->info().volume() = newVol;
- if (EpsVolPercent_RTRG>0) { //if the criterion is not used, then we skip these updates a save a LOT of time when Nthreads > 1
+ if (defTolerance>0) { //if the criterion is not used, then we skip these updates a save a LOT of time when Nthreads > 1
#pragma omp atomic
totVol+=newVol;
#pragma omp atomic
totDVol+=abs(dVol);}
}
- if (EpsVolPercent_RTRG>0) eps_vol_max = totDVol/totVol;
+ if (defTolerance>0) eps_vol_max = totDVol/totVol;
for (unsigned int n=0; n<flow->imposedF.size();n++) {
flow->IFCells[n]->info().dv()+=flow->imposedF[n].second;
flow->IFCells[n]->info().Pcondition=false;}
@@ -508,38 +469,14 @@
template<class Cellhandle>
Real FlowEngine::Volume_cell_single_fictious ( Cellhandle cell )
{
- #if 0
- //Without buffer
- Vector3r V[3];
- int b=0;
- int w=0;
- cell->info().volumeSign=1;
- Real Wall_coordinate=0;
-
- for ( int y=0;y<4;y++ ) {
- if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( y )->info().id(), scene );
- V[w]=sph->state->pos;
- w++;
- } else {
- b = cell->vertex ( y )->info().id();
- const shared_ptr<Body>& wll = Body::byId ( b , scene );
- if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wall_thickness/2;
- else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
- }
- }
- Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
- #else
- //With buffer
- Vector3r V[3];
- int b=0;
- int w=0;
- cell->info().volumeSign=1;
- Real Wall_coordinate=0;
-
- for ( int y=0;y<4;y++ ) {
- if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
-// const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( y )->info().id(), scene );
+ Vector3r V[3];
+ int b=0;
+ int w=0;
+ cell->info().volumeSign=1;
+ Real Wall_coordinate=0;
+
+ for ( int y=0;y<4;y++ ) {
+ if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
V[w]=positionBufferCurrent[cell->vertex ( y )->info().id()].pos;
w++;
} else {
@@ -550,8 +487,6 @@
}
}
Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
- #endif
-
return abs ( Volume );
}
template<class Cellhandle>
@@ -570,16 +505,13 @@
if ( cell->vertex ( g )->info().isFictious ) {
b[j] = cell->vertex ( g )->info().id();
coord[j]=solver->boundary ( b[j] ).coordinate;
-// const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wall_thickness/2;
else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
j++;
} else if ( first_sph ) {
-// const shared_ptr<Body>& sph1 = Body::byId ( cell->vertex ( g )->info().id(), scene );
A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;
first_sph=false;
} else {
-// const shared_ptr<Body>& sph2 = Body::byId ( cell->vertex ( g )->info().id(), scene );
B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;;
}
}
@@ -622,18 +554,10 @@
Real FlowEngine::Volume_cell ( Cellhandle cell )
{
static const Real inv6 = 1/6.;
- #if 0
- //Without buffer
- const Vector3r& p0 = Body::byId ( cell->vertex ( 0 )->info().id(), scene )->state->pos;
- const Vector3r& p1 = Body::byId ( cell->vertex ( 1 )->info().id(), scene )->state->pos;
- const Vector3r& p2 = Body::byId ( cell->vertex ( 2 )->info().id(), scene )->state->pos;
- const Vector3r& p3 = Body::byId ( cell->vertex ( 3 )->info().id(), scene )->state->pos;
- #else
const Vector3r& p0 = positionBufferCurrent[cell->vertex ( 0 )->info().id()].pos;
const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos;
const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos;
const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos;
- #endif
Real volume = inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
return volume;
@@ -642,7 +566,6 @@
void FlowEngine::ComputeViscousForces ( Solver& flow )
{
if (viscousShear || normalLubrication || shearLubrication){
-// flow->ComputeEdgesSurfaces(); //only done in buildTriangulation
if ( Debug ) cout << "Application of viscous forces" << endl;
if ( Debug ) cout << "Number of edges = " << flow.Edge_ids.size() << endl;
for ( unsigned int k=0; k<flow.viscousShearForces.size(); k++ ) flow.viscousShearForces[k]=Vector3r::Zero();
@@ -659,7 +582,6 @@
const int& id2 = flow.Edge_ids[i].second;
int hasFictious= Tes.vertex ( id1 )->info().isFictious + Tes.vertex ( id2 )->info().isFictious;
-// if (hasFictious==2) continue;
if (hasFictious>0 or id1==id2) continue;
const shared_ptr<Body>& sph1 = Body::byId ( id1, scene );
const shared_ptr<Body>& sph2 = Body::byId ( id2, scene );
@@ -714,9 +636,6 @@
visc_f = flow.computeShearLubricationForce(deltaShearV,surfaceDist,i,eps,O1O2,meanRad);
else if (viscousShear)
visc_f = flow.computeViscousShearForce ( deltaShearV, i , Rh);
-
- if ( Debug ) cout << "la force visqueuse entre " << id1 << " et " << id2 << "est " << visc_f << endl;
-
if (viscousShear || shearLubrication){
flow.viscousShearForces[id1]+=visc_f;
@@ -729,9 +648,6 @@
flow.viscousBodyStress[id1] += visc_f * O1C_vect.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
flow.viscousBodyStress[id2] += (-visc_f) * O2C_vect.transpose()/ (4.0/3.0 *3.14* pow(r2,3));}
}
-
-
-
/// Compute the normal lubrication force applied on each particle
if (normalLubrication){
deltaNormV = normal.dot(deltaV);
@@ -743,12 +659,9 @@
if (viscousNormalBodyStress){
flow.lubBodyStress[id1] += lub_f * O1C_vect.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
flow.lubBodyStress[id2] += (-lub_f) *O2C_vect.transpose() / (4.0/3.0 *3.14* pow(r2,3));}
- if ( Debug ) cout << "la force normale entre " << id1 << " et " << id2 << "est " << lub_f << endl;
}
}
-
- if(Debug) cout<<"number of viscousShearForce"<<flow.viscousShearForces.size()<<endl;
}
}
@@ -773,18 +686,18 @@
Build_Triangulation(P_zero,solver);
if (solver->errorCode>0) {LOG_INFO("triangulation error, pausing"); Omega::instance().pause(); return;}
Initialize_volumes(solver); backgroundSolver=solver; backgroundCompleted=true;}
-// if ( first ) {Build_Triangulation ( P_zero ); Update_Triangulation = false; Initialize_volumes();}
+// if ( first ) {Build_Triangulation ( P_zero ); updateTriangulation = false; Initialize_volumes();}
timingDeltas->checkpoint("Triangulating");
UpdateVolumes (solver);
Eps_Vol_Cumulative += eps_vol_max;
- if ( (EpsVolPercent_RTRG>0 && Eps_Vol_Cumulative > EpsVolPercent_RTRG) || retriangulationLastIter>PermuteInterval ) {
- Update_Triangulation = true;
+ if ( (defTolerance>0 && Eps_Vol_Cumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval ) {
+ updateTriangulation = true;
Eps_Vol_Cumulative=0;
retriangulationLastIter=0;
ReTrg++;
} else {
- Update_Triangulation = false;
+ updateTriangulation = false;
retriangulationLastIter++;}
timingDeltas->checkpoint("Update_Volumes");
@@ -816,8 +729,8 @@
///End Compute flow and forces
timingDeltas->checkpoint("Applying Forces");
if (multithread && !first) {
- while (Update_Triangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
- if (Update_Triangulation || ellapsedIter>(0.5*PermuteInterval)) {
+ while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
+ if (updateTriangulation || ellapsedIter>(0.5*meshUpdateInterval)) {
if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
if (fluidBulkModulus>0) solver->Interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
solver=backgroundSolver;
@@ -839,12 +752,12 @@
if (Debug && !backgroundCompleted) cerr<<"still computing solver in the background"<<endl;
ellapsedIter++;}
} else {
- if (Update_Triangulation && !first) {
+ if (updateTriangulation && !first) {
cachedCell= Cell(*(scene->cell));
Build_Triangulation (P_zero, solver);
Initialize_volumes(solver);
ComputeViscousForces(*solver);
- Update_Triangulation = false;}
+ updateTriangulation = false;}
}
first=false;
timingDeltas->checkpoint("Ending");
@@ -995,6 +908,35 @@
//call recursively, since the returned cell could be also a ghost (especially if baseCell is a non-periodic type from the external contour
locateCell ( ch,index,baseIndex,flow,++count );
if ( ch==baseCell ) cerr<<"WTF!!"<<endl;
+ //check consistency
+ bool checkC=false;
+ for (int kk=0; kk<4;kk++) if ((!baseCell->vertex(kk)->info().isGhost) && ((!baseCell->vertex(kk)->info().isFictious))) checkC = true;
+ if (checkC) {
+ bool checkV=true;
+ for (int kk=0; kk<4;kk++) {
+ checkV=false;
+ for (int jj=0; jj<4;jj++)
+ if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
+ if (!checkV) {cerr <<"periodicity is broken"<<endl;
+ for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
+ cerr<<" vs. ";
+ for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
+ cerr<<endl;}
+ }
+ } else {
+// bool checkV=true;
+// for (int kk=0; kk<4;kk++) {
+// checkV=false;
+// for (int jj=0; jj<4;jj++)
+// if (baseCell->vertex(kk)->info().id() == ch->vertex(jj)->info().id()) checkV = true;
+// if (!checkV) {cerr <<"periodicity is broken (that's ok probably)"<<endl;
+// for (int jj=0; jj<4;jj++) cerr<<baseCell->vertex(jj)->info().id()<<" ";
+// cerr<<" vs. ";
+// for (int jj=0; jj<4;jj++) cerr<<ch->vertex(jj)->info().id()<<" ";
+// cerr<<endl;}
+// }
+ }
+
base_info.isGhost=true;
base_info._pression=& ( ch->info().p() );
base_info.index=ch->info().index;
@@ -1129,7 +1071,7 @@
if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0)) flow->Interpolate ( flow->T[!flow->currentTes], Tes );
// if ( !first && (useSolver==0 || fluidBulkModulus>0)) flow->Interpolate ( flow->T[!flow->currentTes], flow->T[flow->currentTes] );
- if ( WaveAction ) flow->ApplySinusoidalPressure ( Tes.Triangulation(), Sinus_Amplitude, Sinus_Average, 30 );
+ if ( WaveAction ) flow->ApplySinusoidalPressure ( Tes.Triangulation(), sineMagnitude, sineAverage, 30 );
if ( viscousShear || normalLubrication || shearLubrication) flow->computeEdgesSurfaces();
if ( Debug ) cout << endl << "end buildTri------" << endl << endl;
}
=== modified file 'pkg/dem/FlowEngine.hpp'
--- pkg/dem/FlowEngine.hpp 2013-05-17 17:19:59 +0000
+++ pkg/dem/FlowEngine.hpp 2013-07-26 18:16:04 +0000
@@ -51,7 +51,7 @@
public :
int retriangulationLastIter;
- enum {wall_left=0, wall_right, wall_bottom, wall_top, wall_back, wall_front};
+ enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax};
Vector3r normal [6];
bool currentTes;
int id_offset;
@@ -118,15 +118,14 @@
void Oedometer_Boundary_Conditions();
void Average_real_cell_velocity();
void saveVtk() {solver->saveVtk();}
- vector<Real> AvFlVelOnSph(unsigned int id_sph) {return solver->Average_Fluid_Velocity_On_Sphere(id_sph);}
+ vector<Real> avFlVelOnSph(unsigned int id_sph) {return solver->Average_Fluid_Velocity_On_Sphere(id_sph);}
- void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; Update_Triangulation=true;}
- void PressureProfile(double wallUpY, double wallDownY) {return solver->MeasurePressureProfile(wallUpY,wallDownY);}
- double MeasurePorePressure(Vector3r pos){return solver->MeasurePorePressure(pos[0], pos[1], pos[2]);}
+// void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; updateTriangulation=true;}
+ void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}
+ double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
TPL int getCell(double posX, double posY, double posZ, Solver& flow){return flow->getCell(posX, posY, posZ);}
- double MeasureAveragedPressure(double posY){return solver->MeasureAveragedPressure(posY);}
- double MeasureTotalAveragedPressure(){return solver->MeasureTotalAveragedPressure();}
-
+ double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
+ double averagePressure(){return solver->averagePressure();}
#ifdef LINSOLV
TPL void exportMatrix(string filename,Solver& flow) {if (useSolver==3) flow->exportMatrix(filename.c_str());
else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
@@ -153,7 +152,6 @@
Real _getCellFlux(unsigned int cond) {return getCellFlux(cond,solver);}
Real _getBoundaryFlux(unsigned int boundary) {return getBoundaryFlux(boundary,solver);}
int _getCell(Vector3r pos) {return getCell(pos[0],pos[1],pos[2],solver);}
-
#ifdef LINSOLV
void _exportMatrix(string filename) {exportMatrix(filename,solver);}
void _exportTriplets(string filename) {exportTriplets(filename,solver);}
@@ -171,72 +169,47 @@
((bool,first,true,,"Controls the initialization/update phases"))
((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
((Real, dt, 0,,"timestep [s]"))
-// ((bool,save_vtk,false,,"Enable/disable vtk files creation for visualization"))
((bool,permeability_map,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
- ((bool, save_mgpost, false,,"Enable/disable mgpost file creation"))
- ((bool, slice_pressures, false, ,"Enable/Disable slice pressure measurement"))
- ((bool, velocity_profile, false, ,"Enable/Disable slice velocity measurement"))
- ((bool, consolidation,false,,"Enable/Disable storing consolidation files"))
((bool, slip_boundary, true,, "Controls friction condition on lateral walls"))
((bool,WaveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
- ((double, Sinus_Amplitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied"))
- ((double, Sinus_Average, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
- ((bool, CachedForces, true,,"Des/Activate the cached forces calculation"))
+ ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))
+ ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
((bool, Debug, false,,"Activate debug messages"))
((double, wall_thickness,0.001,,"Walls thickness"))
- ((double,P_zero,0,,"Initial internal pressure for oedometer test"))
+ ((double,P_zero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
((double,Tolerance,1e-06,,"Gauss-Seidel Tolerance"))
((double,Relax,1.9,,"Gauss-Seidel relaxation"))
- ((bool, Update_Triangulation, 0,,"If true the medium is retriangulated. Can be switched on to force retriangulation after some events (else it will be true periodicaly based on :yref:`FlowEngine::EpsVolPercent_RTRG` and :yref:`FlowEngine::PermuteInterval`"))
- ((int,PermuteInterval,1000,,"Maximum number of timesteps between re-triangulation events. See also :yref:`FlowEngine::EpsVolPercent_RTRG`."))
- ((double, eps_vol_max, 0,,"Maximal absolute volumetric strain computed at each iteration"))
- ((double, EpsVolPercent_RTRG,0.01,,"Percentuage of cumulate eps_vol at which retriangulation of pore space is performed"))
- ((double, porosity, 0,,"Porosity computed at each retriangulation"))
- ((bool,compute_K,false,,"Activates permeability measure within a granular sample"))
- ((bool,meanKStat,false,,"Local permeabilities' correction through an optimized threshold"))
- ((bool,clampKValues,true,,"If true, clamp local permeabilities between min/max*globalK threshold. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))
+ ((bool, updateTriangulation, 0,,"If true the medium is retriangulated. Can be switched on to force retriangulation after some events (else it will be true periodicaly based on :yref:`FlowEngine::defTolerance` and :yref:`FlowEngine::meshUpdateInterval`. Of course, it costs CPU time."))
+ ((int,meshUpdateInterval,1000,,"Maximum number of timesteps between re-triangulation events. See also :yref:`FlowEngine::defTolerance`."))
+ ((double, eps_vol_max, 0,(Attr::readonly),"Maximal absolute volumetric strain computed at each iteration. |yupdate|"))
+ ((double, defTolerance,0.05,,"Cumulated deformation threshold for which retriangulation of pore space is performed. If negative, the triangulation update will occure with a fixed frequency on the basis of :yref:`FlowEngine::meshUpdateInterval`"))
+ ((double, porosity, 0,(Attr::readonly),"Porosity computed at each retriangulation |yupdate|"))
+ ((bool,meanKStat,false,,"report the local permeabilities' correction"))
+ ((bool,clampKValues,true,,"If true, clamp local permeabilities in [minKdivKmean,maxKdivKmean]*globalK. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))
((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
- ((double,permeability_factor,0.0,,"permability multiplier"))
- ((double,viscosity,1.0,,"viscosity of fluid"))
- ((double,stiffness, 10000,,"stiffness modulus"))
- ((Real,loadFactor,1.1,,"Load multiplicator for oedometer test"))
- ((double, K, 0,, "Permeability of the sample"))
+ ((double,permeabilityFactor,0.0,,"permability multiplier"))
+ ((double,viscosity,1.0,,"viscosity of the fluid"))
+ ((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))
((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
-// ((std::string,key,"",,"A string appended at the output files, use it to name simulations."))
- ((double, V_d, 0,,"darcy velocity of fluid in sample"))
- ((bool, Flow_imposed_TOP_Boundary, true,, "if false involve pressure imposed condition"))
- ((bool, Flow_imposed_BOTTOM_Boundary, true,, "if false involve pressure imposed condition"))
- ((bool, Flow_imposed_FRONT_Boundary, true,, "if false involve pressure imposed condition"))
- ((bool, Flow_imposed_BACK_Boundary, true,, "if false involve pressure imposed condition"))
- ((bool, Flow_imposed_LEFT_Boundary, true,, "if false involve pressure imposed condition"))
- ((bool, Flow_imposed_RIGHT_Boundary, true,,"if false involve pressure imposed condition"))
- ((double, Pressure_TOP_Boundary, 0,, "Pressure imposed on top boundary"))
- ((double, Pressure_BOTTOM_Boundary, 0,, "Pressure imposed on bottom boundary"))
- ((double, Pressure_FRONT_Boundary, 0,, "Pressure imposed on front boundary"))
- ((double, Pressure_BACK_Boundary, 0,,"Pressure imposed on back boundary"))
- ((double, Pressure_LEFT_Boundary, 0,, "Pressure imposed on left boundary"))
- ((double, Pressure_RIGHT_Boundary, 0,, "Pressure imposed on right boundary"))
- ((Vector3r, topBoundaryVelocity, Vector3r::Zero(),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
+ ((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))
+ ((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, ymax,3,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, zmin,4,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, zmax,5,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+
+ ((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends."))
+ ((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future)."))
+ //FIXME: to be implemented:
+ ((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
- ((int, wallTopId,3,,"Id of top boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
- ((int, wallBottomId,2,,"Id of bottom boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
- ((int, wallFrontId,5,,"Id of front boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
- ((int, wallBackId,4,,"Id of back boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
- ((int, wallLeftId,0,,"Id of left boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
- ((int, wallRightId,1,,"Id of right boundary (default value is ok if aabbWalls are appended BEFORE spheres.)"))
- ((bool, BOTTOM_Boundary_MaxMin, 1,,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
- ((bool, TOP_Boundary_MaxMin, 1,,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
- ((bool, RIGHT_Boundary_MaxMin, 1,,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
- ((bool, LEFT_Boundary_MaxMin, 1,,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
- ((bool, FRONT_Boundary_MaxMin, 1,,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
- ((bool, BACK_Boundary_MaxMin, 1,,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
- ((bool, areaR2Permeability, 1,,"Use corrected formula for permeabilities calculation in flowboundingsphere (areaR2permeability variable)"))
- ((bool, viscousShear, false,,"Compute viscous shear terms as developped by Donia Marzougui"))
- ((bool, shearLubrication, false,,"Compute shear lubrication force as developped by Brule"))
- ((double, eps, 0.00001,,"minimum distance between particles"))
+ ((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))
+ ((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
+ ((bool, viscousShear, false,,"Compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))
+ ((bool, shearLubrication, false,,"Compute shear lubrication force as developped by Brule (FIXME: ref.) "))
+ ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
((bool, pressureForce, true,,"Compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
-
((bool, normalLubrication, false,,"Compute normal lubrication force as developped by Brule"))
((bool, viscousNormalBodyStress, false,,"Compute normal viscous stress applied on each body"))
((bool, viscousShearBodyStress, false,,"Compute shear viscous stress applied on each body"))
@@ -245,19 +218,17 @@
((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
#endif
- ((Real, allDeprecs, 0,,"transitory variable: point removed attributes to this so that older scripts won't crash."))
,
/*deprec*/
((meanK_opt,clampKValues,"the name changed"))
- ((meanK_correction,allDeprecs,"the name changed"))
,,
timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
- for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero();}
- normal[wall_bottom].y()=normal[wall_left].x()=normal[wall_back].z()=1;
- normal[wall_top].y()=normal[wall_right].x()=normal[wall_front].z()=-1;
+ for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero(); wallIds[i]=i;}
+ normal[wall_ymin].y()=normal[wall_xmin].x()=normal[wall_zmin].z()=1;
+ normal[wall_ymax].y()=normal[wall_xmax].x()=normal[wall_zmax].z()=-1;
solver = shared_ptr<FlowSolver> (new FlowSolver);
first=true;
- Update_Triangulation=false;
+ updateTriangulation=false;
eps_vol_max=Eps_Vol_Cumulative=retriangulationLastIter=0;
ReTrg=1;
backgroundCompleted=true;
@@ -269,27 +240,24 @@
.def("clearImposedPressure",&FlowEngine::_clearImposedPressure,"Clear the list of points with pressure imposed.")
.def("clearImposedFlux",&FlowEngine::_clearImposedFlux,"Clear the list of points with flux imposed.")
.def("getCellFlux",&FlowEngine::_getCellFlux,(python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
- .def("getBoundaryFlux",&FlowEngine::_getBoundaryFlux,(python::arg("boundary")),"Get total flux through boundary defined by its body id.")
+ .def("getBoundaryFlux",&FlowEngine::_getBoundaryFlux,(python::arg("boundary")),"Get total flux through boundary defined by its body id.\n\n.. note::\n\t The flux may be not zero even for no-flow condition. This artifact comes from cells which are incident to two or more boundaries (along the edges of the sample, typically). Such flux evaluation on impermeable boundary is just irrelevant, it does not imply that the boundary condition is not applied properly.")
.def("getConstrictions",&FlowEngine::_getConstrictions,(python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
.def("getConstrictionsFull",&FlowEngine::_getConstrictionsFull,(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("saveVtk",&FlowEngine::saveVtk,"Save pressure field in vtk format.")
- .def("AvFlVelOnSph",&FlowEngine::AvFlVelOnSph,(python::arg("Id_sph")),"Compute a sphere-centered average fluid velocity")
+ .def("avFlVelOnSph",&FlowEngine::avFlVelOnSph,(python::arg("Id_sph")),"Compute a sphere-centered average fluid velocity")
.def("fluidForce",&FlowEngine::_fluidForce,(python::arg("Id_sph")),"Return the fluid force on sphere Id_sph.")
.def("shearLubForce",&FlowEngine::_shearLubForce,(python::arg("Id_sph")),"Return the shear lubrication force on sphere Id_sph.")
.def("shearLubTorque",&FlowEngine::_shearLubTorque,(python::arg("Id_sph")),"Return the shear lubrication torque on sphere Id_sph.")
.def("normalLubForce",&FlowEngine::_normalLubForce,(python::arg("Id_sph")),"Return the normal lubrication force on sphere Id_sph.")
.def("bodyShearLubStress",&FlowEngine::_bodyShearLubStress,(python::arg("Id_sph")),"Return the shear lubrication stress on sphere Id_sph.")
.def("bodyNormalLubStress",&FlowEngine::_bodyNormalLubStress,(python::arg("Id_sph")),"Return the normal lubrication stress on sphere Id_sph.")
- .def("setBoundaryVel",&FlowEngine::setBoundaryVel,(python::arg("vel")),"Change velocity on top boundary.")
- .def("PressureProfile",&FlowEngine::PressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
- .def("MeasurePorePressure",&FlowEngine::MeasurePorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
- .def("MeasureAveragedPressure",&FlowEngine::MeasureAveragedPressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
- .def("MeasureTotalAveragedPressure",&FlowEngine::MeasureTotalAveragedPressure,"Measure averaged pore pressure in the entire volume")
-
+ .def("pressureProfile",&FlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+ .def("getPorePressure",&FlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
+ .def("averageSlicePressure",&FlowEngine::averageSlicePressure,(python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
+ .def("averagePressure",&FlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
.def("updateBCs",&FlowEngine::_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",&FlowEngine::emulateAction,"get scene and run action (may be used to manipulate engine outside the main loop).")
+ .def("emulateAction",&FlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
.def("getCell",&FlowEngine::_getCell,(python::arg("pos")),"get id of the cell containing (X,Y,Z).")
-
#ifdef LINSOLV
.def("exportMatrix",&FlowEngine::_exportMatrix,(python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
.def("exportTriplets",&FlowEngine::_exportTriplets,(python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
@@ -305,7 +273,7 @@
if (!flow) LOG_ERROR("no flow defined yet, run at least one iter");
flow->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) );
//force immediate update of boundary conditions
- Update_Triangulation=true;
+ updateTriangulation=true;
return flow->imposedP.size()-1;
}
@@ -367,12 +335,11 @@
Real _getBoundaryFlux(unsigned int boundary) {return getBoundaryFlux(boundary,solver);}
void _updateBCs() {updateBCs(solver);}
- double MeasurePorePressure(Vector3r pos){return solver->MeasurePorePressure(pos[0], pos[1], pos[2]);}
- double MeasureTotalAveragedPressure(){return solver->MeasureTotalAveragedPressure();}
- void PressureProfile(double wallUpY, double wallDownY) {return solver->MeasurePressureProfile(wallUpY,wallDownY);}
+ double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
+ double averagePressure(){return solver->averagePressure();}
+ void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}
int _getCell(Vector3r pos) {return getCell(pos[0],pos[1],pos[2],solver);}
-
#ifdef LINSOLV
void _exportMatrix(string filename) {exportMatrix(filename,solver);}
void _exportTriplets(string filename) {exportTriplets(filename,solver);}
@@ -394,13 +361,14 @@
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_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine,"An engine to solve flow problem in saturated granular media",
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(PeriodicFlowEngine,FlowEngine,"A variant of :yref:`FlowEngine` implementing periodic boundary conditions. The API is very similar.",
((Real,duplicateThreshold, 0.06,,"distance from cell borders that will triger periodic duplication in the triangulation |yupdate|"))
((Vector3r, gradP, Vector3r::Zero(),,"Macroscopic pressure gradient"))
,,
- wallTopId=wallBottomId=wallFrontId=wallBackId=wallLeftId=wallRightId=-1;
+ wallIds=vector<int>(6,-1);
+// wallTopId=wallBottomId=wallFrontId=wallBackId=wallLeftId=wallRightId=-1;
solver = shared_ptr<FlowSolver> (new FlowSolver);
- Update_Triangulation=false;
+ updateTriangulation=false;
eps_vol_max=Eps_Vol_Cumulative=retriangulationLastIter=0;
ReTrg=1;
first=true;
@@ -417,9 +385,9 @@
.def("saveVtk",&PeriodicFlowEngine::saveVtk,"Save pressure field in vtk format.")
.def("imposePressure",&PeriodicFlowEngine::_imposePressure,(python::arg("pos"),python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
.def("getBoundaryFlux",&PeriodicFlowEngine::_getBoundaryFlux,(python::arg("boundary")),"Get total flux through boundary defined by its body id.")
- .def("MeasurePorePressure",&PeriodicFlowEngine::MeasurePorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
- .def("MeasureTotalAveragedPressure",&PeriodicFlowEngine::MeasureTotalAveragedPressure,"Measure averaged pore pressure in the entire volume")
- .def("PressureProfile",&PeriodicFlowEngine::PressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+ .def("getPorePressure",&PeriodicFlowEngine::getPorePressure,(python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
+ .def("averagePressure",&PeriodicFlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
+ .def("pressureProfile",&PeriodicFlowEngine::pressureProfile,(python::arg("wallUpY"),python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
.def("updateBCs",&PeriodicFlowEngine::_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)")
=== added file 'pkg/dem/LudingPM.cpp'
--- pkg/dem/LudingPM.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/LudingPM.cpp 2013-07-08 08:52:03 +0000
@@ -0,0 +1,177 @@
+#include"LudingPM.hpp"
+#include<yade/core/State.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/core/Scene.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+
+YADE_PLUGIN((LudingMat)(LudingPhys)(Ip2_LudingMat_LudingMat_LudingPhys)(Law2_ScGeom_LudingPhys_Basic));
+
+LudingMat::~LudingMat(){}
+
+LudingPhys::~LudingPhys(){}
+
+void Ip2_LudingMat_LudingMat_LudingPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
+ if(interaction->phys) return;
+
+ LudingMat* mat1 = static_cast<LudingMat*>(b1.get());
+ LudingMat* mat2 = static_cast<LudingMat*>(b2.get());
+
+ const Real k11 = mat1->k1; const Real k12 = mat2->k1;
+ const Real kp1 = mat1->kp; const Real kp2 = mat2->kp;
+ const Real kc1 = mat1->kc; const Real kc2 = mat2->kc;
+ const Real G01 = mat1->G0; const Real G02 = mat2->G0;
+ const Real PhiF1 = mat1->PhiF; const Real PhiF2 = mat2->PhiF;
+
+ LudingPhys* phys = new LudingPhys();
+
+ phys->k1 = this->reduced(k11, k12);
+ phys->kp = this->reduced(kp1, kp2);
+ phys->kc = this->reduced(kc1, kc2);
+ phys->PhiF = this->reduced(PhiF1, PhiF2);
+ phys->k2 = 0.0;
+ phys->G0 = this->reduced(G01, G02);
+
+
+ Real a1 = 0.0;
+ Real a2 = 0.0;
+ Sphere* s1=dynamic_cast<Sphere*>(Body::byId(interaction->getId1())->shape.get());
+ Sphere* s2=dynamic_cast<Sphere*>(Body::byId(interaction->getId2())->shape.get());
+ if (s1 and s2) {
+ a1 = s1->radius;
+ a2 = s2->radius;
+ } else if (s1 and not(s2)) {
+ a1 = s1->radius;
+ } else {
+ a2 = s2->radius;
+ }
+
+ if (phys->k1 >= phys->kp) {
+ throw runtime_error("k1 have to be less as kp!"); // [Luding2008], sentence after equation (6); kp = k2^
+ // [Singh2013], sentence after equation (6)
+ }
+
+
+
+ phys->tangensOfFrictionAngle = std::tan(std::min(mat1->frictionAngle, mat2->frictionAngle));
+
+ phys->shearForce = Vector3r(0,0,0);
+ phys->ThetMax = 0.0;
+ phys->ThetNull = 0.0;
+ phys->ThetPMax = phys->kp/(phys->kp-phys->k1)*phys->PhiF*2*a1*a2/(a1+a2); // [Luding2008], equation (7)
+ // [Singh2013], equation (11)
+
+ interaction->phys = shared_ptr<LudingPhys>(phys);
+}
+
+Real Ip2_LudingMat_LudingMat_LudingPhys::reduced(Real a1, Real a2){
+ Real a = (a1?1/a1:0) + (a2?1/a2:0); a = a?1/a:0;
+ return a;
+}
+
+void Law2_ScGeom_LudingPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
+
+ const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
+ LudingPhys& phys=*static_cast<LudingPhys*>(_phys.get());
+
+ const int id1 = I->getId1();
+ const int id2 = I->getId2();
+
+ const Real Theta = geom.penetrationDepth;
+
+ if (Theta<0) {
+ scene->interactions->requestErase(I);
+ return;
+ };
+
+ const BodyContainer& bodies = *scene->bodies;
+
+
+
+
+ Real forceHys = 0.0;
+
+ if (phys.ThetMax/phys.ThetPMax >= 1.0) { // [Luding2008], equation (8)
+ phys.k2 = phys.kp; // [Singh2013], equation (10)
+ } else {
+ phys.k2 = phys.k1 + (phys.kp - phys.k1)*phys.ThetMax/phys.ThetPMax;
+ }
+
+ if (Theta > phys.ThetMax) {
+ phys.ThetMax = Theta;
+ phys.ThetNull = (1.0 - phys.k1/phys.k2)*phys.ThetMax; // [Luding2008], equation over Fig 1
+ // [Singh2013], equation (8)
+ }
+
+ Real k2DeltaTtmp = phys.k2*(Theta - phys.ThetNull);
+ if ( k2DeltaTtmp >= phys.k1*Theta) { // [Luding2008], equation (6)
+ forceHys = phys.k1*Theta; // [Singh2013], equation (6)
+ } else if (k2DeltaTtmp > -phys.kc*Theta and k2DeltaTtmp < phys.k1*Theta) {
+ forceHys = phys.k2*(Theta-phys.ThetNull);
+ } else if (k2DeltaTtmp<=-phys.kc*Theta) {
+ forceHys = -phys.kc*Theta;
+ }
+
+
+
+ //===================================================================
+ //===================================================================
+ //===================================================================
+ // Copy-paste from ViscoElasticPM
+
+ const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
+ const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
+
+ Vector3r& shearForce = phys.shearForce;
+ if (I->isFresh(scene)) shearForce=Vector3r(0,0,0);
+ const Real& dt = scene->dt;
+ shearForce = geom.rotate(shearForce);
+
+
+ // Handle periodicity.
+
+ const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero();
+ const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero();
+
+
+ const Vector3r c1x = (geom.contactPoint - de1.pos);
+ const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
+
+
+ const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
+ const Real normalVelocity = geom.normal.dot(relativeVelocity);
+ const Vector3r shearVelocity = relativeVelocity-normalVelocity*geom.normal;
+
+ shearForce += phys.ks*dt*shearVelocity; // the elastic shear force have a history, but
+ Vector3r shearForceVisc = Vector3r::Zero(); // the viscous shear damping haven't a history because it is a function of the instant velocity
+
+
+ phys.normalForce = (forceHys + phys.G0 * normalVelocity)*geom.normal;
+
+
+ const Real maxFs = phys.normalForce.squaredNorm() * std::pow(phys.tangensOfFrictionAngle,2);
+ if( shearForce.squaredNorm() > maxFs )
+ {
+ const Real ratio = sqrt(maxFs) / shearForce.norm();
+ shearForce *= ratio;
+ }
+ else
+ {
+ // shearForceVisc = phys.cs*shearVelocity; //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ shearForceVisc = phys.G0*shearVelocity;
+ }
+ //===================================================================
+ //===================================================================
+ //===================================================================
+
+
+ if (I->isActive) {
+ const Vector3r f = phys.normalForce + shearForce + shearForceVisc;
+ addForce (id1,-f,scene);
+ addForce (id2, f,scene);
+ addTorque(id1,-c1x.cross(f),scene);
+ addTorque(id2, c2x.cross(f),scene);
+ }
+
+
+}
=== added file 'pkg/dem/LudingPM.hpp'
--- pkg/dem/LudingPM.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/LudingPM.hpp 2013-07-08 08:52:03 +0000
@@ -0,0 +1,65 @@
+#pragma once
+
+#include<yade/core/Material.hpp>
+#include<yade/pkg/dem/FrictPhys.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+
+
+class LudingMat : public Material {
+ public:
+ virtual ~LudingMat();
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(LudingMat,Material,"Material for simple Ludning`s model of contact.\n",
+ ((Real,k1,NaN,,"Slope of loading plastic branch"))
+ ((Real,kp,NaN,,"Slope of unloading and reloading limit elastic branch"))
+ ((Real,kc,NaN,,"Slope of irreversible, tensile adhesive branch"))
+ ((Real,PhiF,NaN,,"Dimensionless plasticity depth"))
+ ((Real,G0,NaN,,"Viscous damping"))
+ ((Real,frictionAngle,NaN,,"Friction angle [rad]")),
+ createIndex();
+ );
+ REGISTER_CLASS_INDEX(LudingMat,Material);
+};
+REGISTER_SERIALIZABLE(LudingMat);
+
+class LudingPhys : public FrictPhys{
+ public:
+ virtual ~LudingPhys();
+ Real R;
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(LudingPhys,FrictPhys,"IPhys created from :yref:`LudingMat`, for use with :yref:`Law2_ScGeom_LudingPhys_Basic`.",
+ ((Real,k1,NaN,,"Slope of loading plastic branch"))
+ ((Real,k2,NaN,,"Slope of unloading and reloading elastic branch"))
+ ((Real,kp,NaN,,"Slope of unloading and reloading limit elastic branch"))
+ ((Real,kc,NaN,,"Slope of irreversible, tensile adhesive branch"))
+ ((Real,PhiF,NaN,,"Dimensionless plasticity depth"))
+ ((Real,ThetMax,NaN,,"Maximum overlap between particles for a collision"))
+ ((Real,ThetPMax,NaN,,"Maximum overlap between particles for the limit case"))
+ ((Real,ThetNull,NaN,,"Force free overlap, plastic contact deformation"))
+ ((Real,G0,NaN,,"Viscous damping")),
+ createIndex();
+ )
+};
+REGISTER_SERIALIZABLE(LudingPhys);
+
+class Ip2_LudingMat_LudingMat_LudingPhys: public IPhysFunctor {
+ public :
+ virtual void go(const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ YADE_CLASS_BASE_DOC(Ip2_LudingMat_LudingMat_LudingPhys,IPhysFunctor,"Convert 2 instances of :yref:`LudingMat` to :yref:`LudingPhys` using the rule of consecutive connection.");
+ FUNCTOR2D(LudingMat,LudingMat);
+ private:
+ Real reduced(Real, Real);
+
+};
+REGISTER_SERIALIZABLE(Ip2_LudingMat_LudingMat_LudingPhys);
+
+class Law2_ScGeom_LudingPhys_Basic: public LawFunctor {
+ public :
+ virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ private:
+ Real calculateCapillarForce(const ScGeom& geom, LudingPhys& phys);
+ FUNCTOR2D(ScGeom,LudingPhys);
+ YADE_CLASS_BASE_DOC(Law2_ScGeom_LudingPhys_Basic,LawFunctor,"Linear viscoelastic model operating on :yref:`ScGeom` and :yref:`LudingPhys`.");
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_LudingPhys_Basic);
=== removed file 'pkg/dem/ParticleSizeDistrbutionRPMRecorder.cpp'
--- pkg/dem/ParticleSizeDistrbutionRPMRecorder.cpp 2010-10-18 14:51:14 +0000
+++ pkg/dem/ParticleSizeDistrbutionRPMRecorder.cpp 1970-01-01 00:00:00 +0000
@@ -1,293 +0,0 @@
-#include"ParticleSizeDistrbutionRPMRecorder.hpp"
-
-YADE_PLUGIN((ParticleSizeDistrbutionRPMRecorder));
-CREATE_LOGGER(ParticleSizeDistrbutionRPMRecorder);
-
-void ParticleSizeDistrbutionRPMRecorder::action() {
- numberCohesiveContacts=0;
- int curBin = 1; //Current bin number for PSD
- vector<identicalIds> arrayIdentIds;
- arrayIdentIds.clear();
-
- FOREACH(const shared_ptr<Body>& b, *scene->bodies){ // Annulling specimenNumber
- if (!b) continue;
- YADE_PTR_CAST<RpmState>(b->state)->specimenNumber = 0;
- YADE_PTR_CAST<RpmState>(b->state)->specimenMass = 0;
- YADE_PTR_CAST<RpmState>(b->state)->specimenVol = 0;
- YADE_PTR_CAST<RpmState>(b->state)->specimenMaxDiam = 0;
- }
-
- //Check all interactions
- FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){
- if(!i->isReal()) continue; //Check whether they are real
- const shared_ptr<RpmPhys>& contPhys = YADE_PTR_CAST<RpmPhys>(i->phys);
-
- Body::id_t id1 = i->getId1(); //Get bodies ids from interaction
- Body::id_t id2 = i->getId2();
-
- const Sphere* sphere1 = dynamic_cast<Sphere*>(Body::byId(id1)->shape.get());
- const Sphere* sphere2 = dynamic_cast<Sphere*>(Body::byId(id2)->shape.get());
-
- if ((sphere1)&&(sphere2)) { //This class is ONLY for spheres
- int& specimenNumberId1 = YADE_PTR_CAST<RpmState>(Body::byId(id1)->state)->specimenNumber;
- int& specimenNumberId2 = YADE_PTR_CAST<RpmState>(Body::byId(id2)->state)->specimenNumber;
-
- bool cohesiveState = contPhys->isCohesive;
- if (cohesiveState==true){
- if ((specimenNumberId1 == 0) and (specimenNumberId2 == 0)){ //Both bodies are "untouched"
- specimenNumberId1 = curBin;
- specimenNumberId2 = curBin;
- curBin++;
- } else if ((specimenNumberId1 != 0) and (specimenNumberId2 == 0)){ //On of bodies is 0, another one has number specimen
- specimenNumberId2 = specimenNumberId1;
- } else if ((specimenNumberId1 == 0) and (specimenNumberId2 != 0)){ //On of bodies is 0, another one has number specimen
- specimenNumberId1 = specimenNumberId2;
- } else if ((specimenNumberId1 != 0) and (specimenNumberId2 != 0) and (specimenNumberId1 != specimenNumberId2) ){ //Bodies have different specimen number, but they are cohesive! Update it
- int minIdR = std::min(specimenNumberId1, specimenNumberId2);
- int maxIdR = std::max(specimenNumberId1, specimenNumberId2);
- specimenNumberId1 = minIdR;
- specimenNumberId2 = minIdR;
- identicalIds tempVar(minIdR, maxIdR, 0, 0);
- arrayIdentIds.push_back(tempVar); //Put in the container, that 2 ids belong to the same spcimen!
- }
- }
- if (contPhys->isCohesive==true) { //Check whether they are cohesive
- numberCohesiveContacts++; //If yes - calculate them
- }
- }
- }
-
- //Clean dublicates
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- for (unsigned int d=0; d<arrayIdentIds.size(); d++) {
- if ((i!=d)&&(arrayIdentIds[i].id1 == arrayIdentIds[d].id1)&&(arrayIdentIds[i].id2 == arrayIdentIds[d].id2)) {
- arrayIdentIds.erase(arrayIdentIds.begin()+d);
- }
- }
- }
-
- //Updating the container.
- bool flagChange=true;
- while (flagChange){
- flagChange=false;
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- for (unsigned int d=0; d<arrayIdentIds.size(); d++) {
- if (i!=d){
- identicalIds& tempAr1 = arrayIdentIds[i];
- identicalIds& tempAr2 = arrayIdentIds[d]; // Initial New
- if ((tempAr1.id2 == tempAr2.id1) and (tempAr2.id1>tempAr1.id1)) { // 1 = 2 |=> 1 = 2
- tempAr2.id1 = tempAr1.id1; // 2 = 3 |=> 1 = 3
- flagChange=true;
- }
- if (tempAr1.id2 == tempAr2.id2) { // Initial New
- if (tempAr2.id1>tempAr1.id1) { // 1 = 3 |=> 1 = 3
- tempAr2.id2 = tempAr2.id1; // 2 = 3 |=> 1 = 2
- tempAr2.id1 = tempAr1.id1;
- flagChange=true;
- } else if (tempAr2.id1<tempAr1.id1) { // Initial New
- tempAr1.id2 = tempAr1.id1; // 2 = 3 |=> 1 = 2
- tempAr1.id1 = tempAr2.id1; // 1 = 3 |=> 1 = 3
- flagChange=true;
- } else { // Initial New
- arrayIdentIds.erase(arrayIdentIds.begin()+d); // 1 = 3 |=> 1 = 3
- } // 1 = 3 |=> DELETE
- }
- }
- }
- }
- }
-
- //Update RpmState, specimenIds
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if (!b) continue;
- const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
- if (sphere) {
- int tempSpecimenNum = YADE_PTR_CAST<RpmState>(b->state)->specimenNumber;
- if ((tempSpecimenNum==arrayIdentIds[i].id1) or (tempSpecimenNum==arrayIdentIds[i].id2)){
- YADE_PTR_CAST<RpmState>(b->state)->specimenNumber=arrayIdentIds[i].id1;
- }
- }
- }
- }
-
- int maximalSpecimenId = curBin;
- arrayIdentIds.clear();
- Real totalMass = 0;
- Real totalVol = 0;
- Real totalPartNum = 0;
- Real const constForVol = 4.0/3.0;
- //Calculate specimen masses, create vector for storing it
- FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if (!b) continue;
- const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
- if (sphere) {
- Real massTemp = b->state->mass;
- Real volTemp = constForVol*Mathr::PI*pow ( sphere->radius, 3 );
- totalMass += massTemp;
- totalVol += volTemp;
- totalPartNum += 1;
- int specimenNumberId = YADE_PTR_CAST<RpmState>(b->state)->specimenNumber;
-
- if (specimenNumberId != 0) { //Check, whether particle already belongs to any specimen
- bool foundItemInArray = false; //If yes, search for suitable "bin"
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- if (arrayIdentIds[i].id1 == specimenNumberId) {
- arrayIdentIds[i].mass+=massTemp;//If "bin" for particle with this specimenId found, put its mass there
- arrayIdentIds[i].vol+=volTemp;//If "bin" for particle with this specimenId found, put its volume there
- arrayIdentIds[i].particleNumber+=1; //Calculate the number of particles in on specimen
- foundItemInArray = true;
- }
- if (foundItemInArray) break;
- }
- if (!foundItemInArray) { //If "bin" for particle is not found, create a "bin" for it
- identicalIds tempVar(specimenNumberId, specimenNumberId+1, massTemp, volTemp);
- arrayIdentIds.push_back(tempVar);
- }
- } else { //If the particle was not in contact with other bodies, we give it maximal possible specimenId
- YADE_PTR_CAST<RpmState>(b->state)->specimenNumber = maximalSpecimenId;
- identicalIds tempVar(maximalSpecimenId, maximalSpecimenId+1, massTemp, volTemp);
- arrayIdentIds.push_back(tempVar);
- maximalSpecimenId++;
- }
- }
- }
-
-
- //Find maximal distance between spheres of one specimen
- FOREACH(const shared_ptr<Body>& b1, *scene->bodies){
- if (!b1) continue;
- FOREACH(const shared_ptr<Body>& b2, *scene->bodies){
- if (!b2) continue;
- const Sphere* sphere1 = dynamic_cast<Sphere*>(b1->shape.get()); //Check, whether it is a sphere
- const Sphere* sphere2 = dynamic_cast<Sphere*>(b2->shape.get());
- int specimenNumberId1 = YADE_PTR_CAST<RpmState>(b1->state)->specimenNumber; //Get specimenNumberId
- int specimenNumberId2 = YADE_PTR_CAST<RpmState>(b2->state)->specimenNumber;
-
- if (((sphere1)&&(sphere2))&&(b1 != b2)&&(specimenNumberId1==specimenNumberId2)) {
- Real distBetweenSpheres = (b1->state->pos - b2->state->pos).norm() + sphere1->radius + sphere2->radius;
- Real maxX = abs(b1->state->pos[0] - b2->state->pos[0]) + sphere1->radius + sphere2->radius;
- Real maxY = abs(b1->state->pos[1] - b2->state->pos[1]) + sphere1->radius + sphere2->radius;
- Real maxZ = abs(b1->state->pos[2] - b2->state->pos[2]) + sphere1->radius + sphere2->radius;
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- if ((arrayIdentIds[i].id1 == specimenNumberId1) or (arrayIdentIds[i].id1 == specimenNumberId2)) {
- if (arrayIdentIds[i].maxDistanceBetweenSpheres<distBetweenSpheres) {
- arrayIdentIds[i].maxDistanceBetweenSpheres = distBetweenSpheres;
- }
- if (arrayIdentIds[i].maxX < maxX) {arrayIdentIds[i].maxX = maxX;}
- if (arrayIdentIds[i].maxY < maxY) {arrayIdentIds[i].maxY = maxY;}
- if (arrayIdentIds[i].maxZ < maxZ) {arrayIdentIds[i].maxZ = maxZ;}
- break;
- }
- }
- }
- }
- }
-
- //Update specimen masses and volumes
- FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if (!b) continue;
- const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
- if (sphere) {
- int specimenNumberId = YADE_PTR_CAST<RpmState>(b->state)->specimenNumber;
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- if (arrayIdentIds[i].id1 == specimenNumberId) {
- YADE_PTR_CAST<RpmState>(b->state)->specimenMass = arrayIdentIds[i].mass; //Each particle will contain now the mass of specimen, to which it belongs to
- if (arrayIdentIds[i].maxDistanceBetweenSpheres==0) {
- arrayIdentIds[i].maxDistanceBetweenSpheres=sphere->radius*2.0;
-
- arrayIdentIds[i].maxX = sphere->radius*2.0;
- arrayIdentIds[i].maxY = sphere->radius*2.0;
- arrayIdentIds[i].maxZ = sphere->radius*2.0;
- }
- YADE_PTR_CAST<RpmState>(b->state)->specimenMaxDiam = arrayIdentIds[i].maxDistanceBetweenSpheres; //Each particle will contain now the maximal diametr of the specimen, to which it belongs to
- YADE_PTR_CAST<RpmState>(b->state)->specimenVol = arrayIdentIds[i].vol; //Each particle will contain now the volume of the specimen, to which it belongs to
- break;
- }
- }
- }
- }
- std::sort (arrayIdentIds.begin(), arrayIdentIds.end(), identicalIds::sortArrayIdentIds);
-
- //Material Analyze===============================================================================================
- vector<materialAnalyze> materialAnalyzeIds;
- materialAnalyzeIds.clear();
- FOREACH(const shared_ptr<Body>& b, *scene->bodies){
- if (!b) continue;
- const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
- if (sphere) {
- int specimenNumberId = YADE_PTR_CAST<RpmState>(b->state)->specimenNumber;
- int materialId = b->material->id;
-
- //Check, whether this specimenId is in array
- bool foundSuitableRecording = false;
- for (unsigned int i=0; i<materialAnalyzeIds.size(); i++) {
- if ((materialAnalyzeIds[i].specId == specimenNumberId) and (materialAnalyzeIds[i].matId == materialId)) {
- materialAnalyzeIds[i].particleNumber++;
- materialAnalyzeIds[i].mass+=b->state->mass;
- materialAnalyzeIds[i].vol+=constForVol*Mathr::PI*pow ( sphere->radius, 3 );
- foundSuitableRecording = true;
- break;
- }
- }
- //If not found the recording, create one
- if (!foundSuitableRecording) {
- Real spheresVolume = constForVol*Mathr::PI*pow ( sphere->radius, 3 );
- materialAnalyze tempVar (materialId, specimenNumberId, 1, b->state->mass, spheresVolume);
- materialAnalyzeIds.push_back(tempVar);
- }
- }
- }
- std::sort (materialAnalyzeIds.begin(), materialAnalyzeIds.end(), materialAnalyze::sortMaterialAnalyze);
-
- //Define, how many material columns we need:
- vector<int> materialCount;
- materialCount.clear();
- for (unsigned int i=0; i<materialAnalyzeIds.size(); i++) {
- bool foundItem = false;
- for (unsigned int w=0; w<materialCount.size(); w++) {
- if (materialCount[w]==materialAnalyzeIds[i].matId) {
- foundItem = true;
- break;
- }
- }
- if (foundItem==false) {materialCount.push_back(materialAnalyzeIds[i].matId);}
- }
- std::sort (materialCount.begin(), materialCount.end());
-
- //=================================================================================================================
- //Save data to a file
- out<<"**********\n";
- out<<"iter\ttotalMass\ttotalVol\ttotalPartNum\tnumSpec\tmatNum\n";
- out<<scene->iter<<"\t"<<totalMass<<"\t"<<totalVol<<"\t"<<totalPartNum<<"\t"<<arrayIdentIds.size()<<"\t"<<materialCount.size()<<"\n\n";
- out<<"id\tmassSpec\tvolSpec\tmaxDiamSpec\tmaxX\tmaxY\tmaxZ\tpartNum\t";
-
-
- if (materialCount.size() > 1) {
- for (unsigned int w=0; w<materialCount.size(); w++) { out<<"mat_"<<materialCount[w]<<"_Mass\tmat_"<<materialCount[w]<<"_Vol\tmat_"<<materialCount[w]<<"_PartNum\t";}
- }
- out<<"\n";
-
- for (unsigned int i=0; i<arrayIdentIds.size(); i++) {
- out<<arrayIdentIds[i].id1<<"\t"<<arrayIdentIds[i].mass<<"\t"<<arrayIdentIds[i].vol<<"\t"<<arrayIdentIds[i].maxDistanceBetweenSpheres<<"\t"<<arrayIdentIds[i].maxX<<"\t"<<arrayIdentIds[i].maxY<<"\t"<<arrayIdentIds[i].maxZ<<"\t"<<arrayIdentIds[i].particleNumber<<"\t";
- if (materialCount.size() > 1) {
- //Find Material Info
- for (unsigned int w=0; w<materialCount.size(); w++) {
- bool findItem=false;
- vector<int> indexToDelete;
- for (unsigned int l=0; l<materialAnalyzeIds.size(); l++) {
- if ((materialAnalyzeIds[l].matId==materialCount[w])&&(materialAnalyzeIds[l].specId==arrayIdentIds[i].id1)) {
- out<<materialAnalyzeIds[l].mass<<"\t"<<materialAnalyzeIds[l].vol<<"\t"<<materialAnalyzeIds[l].particleNumber<<"\t";
- findItem=true;
- }
- }
- if (!findItem) {
- out<<"0\t0\t0\t";
- }
- }
- }
- out<<"\n";
- }
-
- out.close();
-}
=== removed file 'pkg/dem/ParticleSizeDistrbutionRPMRecorder.hpp'
--- pkg/dem/ParticleSizeDistrbutionRPMRecorder.hpp 2010-11-12 08:03:16 +0000
+++ pkg/dem/ParticleSizeDistrbutionRPMRecorder.hpp 1970-01-01 00:00:00 +0000
@@ -1,51 +0,0 @@
-/*! This class is for recording Particle Size Distribution
- * of RPM model in a file. Class derived from Recorder
-*/
-#pragma once
-#include <yade/pkg/common/Recorder.hpp>
-#include <yade/pkg/common/Dispatching.hpp>
-#include <yade/pkg/dem/RockPM.hpp>
-#include<yade/core/Scene.hpp>
-#include<yade/pkg/common/Sphere.hpp>
-
-class ParticleSizeDistrbutionRPMRecorder: public Recorder {
- std::ofstream outFile;
- public:
- virtual void action();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(ParticleSizeDistrbutionRPMRecorder,Recorder,
- "Store number of PSD in RPM model to file.",
- ((int,numberCohesiveContacts,0,,"Number of cohesive contacts found at last run. [-]")),
- initRun=true;);
- DECLARE_LOGGER;
-};
-
-struct identicalIds{
- int id1, id2,particleNumber;
- Real mass, vol, maxDistanceBetweenSpheres, maxX, maxY, maxZ;
- identicalIds (int id1r, int id2r, Real massr, Real volr){
- assert(id1r<id2r);
- id1 = id1r;
- id2 = id2r;
- mass = massr;
- vol = volr;
- maxDistanceBetweenSpheres = 0;
- particleNumber = 1;
- maxX = 0; maxY = 0; maxZ = 0;
- }
- static bool sortArrayIdentIds (identicalIds i, identicalIds d) {return i.mass>d.mass;}
-};
-
-struct materialAnalyze{
- int matId, specId, particleNumber;
- Real mass, vol;
- materialAnalyze (int matIdR, int specIdR, int particleNumberR, Real massR, Real volR){
- matId = matIdR;
- specId = specIdR;
- particleNumber = particleNumberR;
- mass = massR;
- vol = volR;
- }
- static bool sortMaterialAnalyze (materialAnalyze i, materialAnalyze d) {return d.specId>i.specId;}
-};
-
-REGISTER_SERIALIZABLE(ParticleSizeDistrbutionRPMRecorder);
=== modified file 'pkg/dem/PeriIsoCompressor.cpp'
--- pkg/dem/PeriIsoCompressor.cpp 2013-05-03 16:45:52 +0000
+++ pkg/dem/PeriIsoCompressor.cpp 2013-07-09 12:42:00 +0000
@@ -264,15 +264,16 @@
}
// variables used in evaluation of ideal stress and ideal strain for each part defined by ##Path
- paths[0]=&xxPath; paths[1]=&yyPath; paths[2]=&zzPath; paths[3]=&yzPath; paths[4]=&zxPath; paths[5]=&xyPath; // pointers to the Paths
+ //paths[0]=&xxPath; paths[1]=&yyPath; paths[2]=&zzPath; paths[3]=&yzPath; paths[4]=&zxPath; paths[5]=&xyPath; // pointers to the Paths
pathSizes[0]=xxPath.size(); pathSizes[1]=yyPath.size(); pathSizes[2]=zzPath.size();
pathSizes[3]=yzPath.size(); pathSizes[4]=zxPath.size(); pathSizes[5]=xyPath.size();
for (int i=0; i<6; i++) {pathsCounter[i] = 0;} // inidicator in which part of the path we are
- // path[0] is a pointer to xxPath
- // path[0]->operator[](j) is j-th Vector2r in path[0]
+ // abcPath[j] is j-th Vector2r in path[0]
// PATH_OP_OP(0,j,k) = path[0]->operator[](j).operator(k) is k-th element of j-th Vector2r of xxPath
- #define PATH_OP_OP(pi,op1i,op2i) paths[pi]->operator[](op1i).operator()(op2i)
+ //#define PATH_OP_OP(pi,op1i,op2i) paths[pi]->operator[](op1i).operator()(op2i)
+ //#define PATH_OP_OP(pi,op1i,op2i) (pi==0?xxPath:pi==1?yyPath:pi==2?zzPath:pi==3?yzPath:pi==4?zxPath:pi==5?xyPath:NULL)->operator[](op1i).operator()(op2i)
+ #define PATH_OP_OP(pi,op1i,op2i) (pi==0?xxPath:pi==1?yyPath:pi==2?zzPath:pi==3?yzPath:pi==4?zxPath:xyPath)[op1i].operator()(op2i)
for (int i=0; i<6; i++) {
for (int j=1; j<pathSizes[i]; j++) {
=== modified file 'pkg/dem/PeriIsoCompressor.hpp'
--- pkg/dem/PeriIsoCompressor.hpp 2012-09-08 01:19:45 +0000
+++ pkg/dem/PeriIsoCompressor.hpp 2013-07-09 12:49:41 +0000
@@ -79,14 +79,14 @@
class Peri3dController: public BoundaryController{
public:
- Vector6r stressOld, stressGoal, strainGoal;
- Vector6i pe, ps;
+ Vector6r stressOld;//, stressGoal, strainGoal;
+ //Vector6i pe, ps;
Matrix3r sigma, epsilon, epsilonRate, rot, nonrot;
- int pathSizes[6], pathsCounter[6], lenPe, lenPs;
- vector<Vector2r>* paths[6];
+ //int pathSizes[6], pathsCounter[6], lenPe, lenPs;
+ //vector<Vector2r>* paths[6];
virtual void action();
- YADE_CLASS_BASE_DOC_ATTRS(Peri3dController,BoundaryController,"Experimental controller of full strain/stress tensors on periodic cell. Detailed documentation is in py/_extraDocs.py.",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Peri3dController,BoundaryController,"Experimental controller of full strain/stress tensors on periodic cell. Detailed documentation is in py/_extraDocs.py.",
((Vector6r,stress,Vector6r::Zero(),,"Current stress vector ($\\sigma_x$,$\\sigma_y$,$\\sigma_z$,$\\tau_{yz}$,$\\tau_{zx}$,$\\tau_{xy}$)|yupdate|."))
((Vector6r,strain,Vector6r::Zero(),,"Current strain (deformation) vector ($\\varepsilon_x$,$\\varepsilon_y$,$\\varepsilon_z$,$\\gamma_{yz}$,$\\gamma_{zx}$,$\\gamma_{xy}$) |yupdate|."))
((Vector6r,strainRate,Vector6r::Zero(),,"Current strain rate vector."))
@@ -108,9 +108,22 @@
((Real,maxStrain,1e6,,"Maximal asolute value of :yref:`strain<Peri3dController.strain>` allowed in the simulation. If reached, the simulation is considered as finished"))
((Real,youngEstimation,1e20,,"Estimation of macroscopic Young's modulus, used for the first simulation step"))
((Real,poissonEstimation,.25,,"Estimation of macroscopic Poisson's ratio, used used for the first simulation step"))
+ //
+ ((Vector6r,stressGoal,Vector6r::Zero(),Attr::readonly,"Peri3dController internal variable"))
+ ((Vector6r,strainGoal,Vector6r::Zero(),Attr::readonly,"Peri3dController internal variable"))
+ ((Vector6i,pe,Vector6i::Zero(),Attr::readonly,"Peri3dController internal variable"))
+ ((Vector6i,ps,Vector6i::Zero(),Attr::readonly,"Peri3dController internal variable"))
+ ((Vector6i,pathSizes,Vector6i::Zero(),Attr::readonly,"Peri3dController internal variable"))
+ ((Vector6i,pathsCounter,Vector6i::Zero(),Attr::readonly,"Peri3dController internal variable"))
+ ((int,lenPe,NaN,Attr::readonly,"Peri3dController internal variable"))
+ ((int,lenPs,NaN,Attr::readonly,"Peri3dController internal variable"))
// not yet used
//((Real,currUnbalanced,NaN,,"current unbalanced force |yupdate|"))
//((Real,maxUnbalanced,1e-4,,"Maximum unbalanced force"))
+ ,
+ /*ctor*/
+ ,
+ /*py*/
);
DECLARE_LOGGER;
=== removed file 'pkg/dem/RockPM.cpp'
--- pkg/dem/RockPM.cpp 2012-09-14 20:18:27 +0000
+++ pkg/dem/RockPM.cpp 1970-01-01 00:00:00 +0000
@@ -1,202 +0,0 @@
-#include"RockPM.hpp"
-#include<yade/core/Scene.hpp>
-#include<yade/pkg/dem/DemXDofGeom.hpp>
-#include<yade/pkg/dem/Shop.hpp>
-#include<math.h>
-
-YADE_PLUGIN((RpmState)(Law2_Dem3DofGeom_RockPMPhys_Rpm)(RpmMat)(Ip2_RpmMat_RpmMat_RpmPhys)(RpmPhys));
-
-
-/********************** Law2_Dem3DofGeom_RockPMPhys_Rpm ****************************/
-CREATE_LOGGER(Law2_Dem3DofGeom_RockPMPhys_Rpm);
-
-void Law2_Dem3DofGeom_RockPMPhys_Rpm::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
- Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
- RpmPhys* phys=static_cast<RpmPhys*>(ip.get());
- const BodyContainer& bodies = *scene->bodies;
-
- const Real& crossSection=phys->crossSection;
-
- Real& epsN=phys->epsN;
- Vector3r& epsT=phys->epsT;
-
- Real sigmaN=phys->sigmaN;
- Vector3r sigmaT=phys->sigmaT;
-
- epsN=geom->strainN();
- epsT=geom->strainT();
-
- const Real& E = phys->E;
- const Real& G = phys->G;
-
- Real& Fn=phys->Fn;
- Vector3r& Fs=phys->Fs;
-
- sigmaN=E*epsN;
- sigmaT=G*epsT;
-
- const int id1 = contact->getId1();
- const int id2 = contact->getId2();
-
- const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
- const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
-
- const Vector3r c1x = (geom->contactPoint - de1.pos);
- const Vector3r c2x = (geom->contactPoint - de2.pos);
-
- const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) ;
- const Real normalVelocity = geom->normal.dot(relativeVelocity);
- const Vector3r shearVelocity = relativeVelocity-normalVelocity*geom->normal;
-
- if (phys->isCohesive) {
- if (epsN<-phys->epsMaxCompression) {
- //std::cout<<"Destruction from Compression\n";
- phys->isCohesive = false;
- } else if (epsN>phys->epsMaxTension) {
- //std::cout<<"Destruction from Tension\n";
- phys->isCohesive = false;
- } else if (epsT.squaredNorm()>(phys->epsMaxShear*phys->epsMaxShear)) {
- //std::cout<<"Destruction from Shear\n";
- phys->isCohesive = false;
- }
- }
-
-
- if (phys->isCohesive) { //Cohesive State
- Fn=sigmaN*crossSection; phys->normalForce=Fn*geom->normal-phys->cn*normalVelocity*geom->normal;
- Fs=sigmaT*crossSection; phys->shearForce = Fs-phys->cs*shearVelocity;
- } else { //UnCohesive State
- if (epsN<0) { //Compression
- Fn=sigmaN*crossSection;
- phys->normalForce=Fn*geom->normal-phys->cn*normalVelocity*geom->normal; //Normal Forces
-
- Fs=sigmaT*crossSection;
- Real maxFs = phys->normalForce.norm()*phys->tanFrictionAngle; /**Check for maximal possible friction force**/
- if(Fs.norm()>maxFs) { Fs*=maxFs/Fs.norm();}
-
- phys->shearForce = Fs-phys->cs*shearVelocity;
-
- } else if (epsN>0) { //If UnCohesive State and particles are too far from each other, delete the interaction
- scene->interactions->requestErase(id1,id2);
- return;
- }
- }
-
-
- //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- //~~~~~~~~~~Destruction Experimental~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- #if 0
- //EXPERIMENTAL!!!
- if (phys->isCohesive) {
- if ((epsN<0) && (-epsN/phys->epsMaxCompression + epsT.norm()/phys->epsMaxShear) > 1) { phys->isCohesive = false;}
- else if ((epsN>0) && (epsN/phys->epsMaxTension + epsT.norm()/phys->epsMaxShear) > 1) { phys->isCohesive = false;}
- }
- #endif
- //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- //~~~~~~~~~~Destruction Experimental~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- #if 0
- sigmaT=G*epsT;
- Real maxFsSq = phys->normalForce.squaredNorm()*phys->tanFrictionAngle*phys->tanFrictionAngle; /**Check for maximal possible friction force**/
- if(sigmaT.squaredNorm()>maxFsSq) {
- sigmaT*=sqrt(maxFsSq/(sigmaT.squaredNorm()));
- }
- Fs=sigmaT*crossSection;
- phys->shearForce = Fs;
- #endif
- //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-
- applyForceAtContactPoint(phys->normalForce + phys->shearForce, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position);
-}
-
-
-CREATE_LOGGER(Ip2_RpmMat_RpmMat_RpmPhys);
-
-void Ip2_RpmMat_RpmMat_RpmPhys::go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction){
- if(interaction->phys) return;
-
- Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->geom.get());
- assert(contGeom);
- //LOG_WARN(Omega::instance().getCurrentIteration());
- const shared_ptr<RpmMat>& rpm1=YADE_PTR_CAST<RpmMat>(pp1);
- const shared_ptr<RpmMat>& rpm2=YADE_PTR_CAST<RpmMat>(pp2);
-
- const int id1 = interaction->getId1();
- const int id2 = interaction->getId2();
-
- Body* b1=Body::byId(id1,scene).get();
- Body* b2=Body::byId(id2,scene).get();
-
- long cohesiveThresholdIter=2;
-
- bool initCohesive = rpm1->initCohesive*rpm2->initCohesive;
-
- Real& E1 = rpm1->young;
- Real& E2 = rpm2->young;
- Real& V1 = rpm1->poisson;
- Real& V2 = rpm2->poisson;
- Real& R1 = contGeom->refR1;
- Real& R2 = contGeom->refR2;
-
- Real E12=2*E1*E2/(E1+E2);
- Real minRad=(R1<=0?R2:(R2<=0?R1:min(R1,R2)));
- Real S12=Mathr::PI*pow(minRad,2);
- Real G_over_E = (rpm1->G_over_E + rpm2->G_over_E)/2;
-
-
- shared_ptr<RpmPhys> contPhys(new RpmPhys());
- contPhys->E=E12;
- contPhys->G=E12*G_over_E;
-
- if ((rpm1->frictionAngle==0) or (rpm2->frictionAngle==0)) { //If one of frictionAngle==0, we have no friction
- contPhys->tanFrictionAngle = 0;
- } else {
- contPhys->tanFrictionAngle=tan(.5*(rpm1->frictionAngle+rpm2->frictionAngle));
- }
-
- contPhys->crossSection=S12;
-
- //We take the minimal strength parameter
- contPhys->epsMaxCompression=(std::min(rpm1->stressCompressMax,rpm2->stressCompressMax))/contPhys->E;
- contPhys->epsMaxTension=(std::min(rpm1->stressStretchingMax,rpm2->stressStretchingMax))/contPhys->E;
- contPhys->epsMaxShear=(std::min(rpm1->stressShearMax,rpm2->stressShearMax))/contPhys->G;
-
- Real RR1 = R1; Real RR2 = R2;
-
- if ((RR1=-1)&&(RR2!=-1)) {
- RR1 = RR2;
- } else if ((RR1!=-1)&&(RR2!=-1)) {
- RR2 = RR1;
- }
- contPhys->Kn = 2*E1*RR1*E2*RR2/(E1*RR1+E2*RR2);
- contPhys->Ks = 2*E1*RR1*V1*E2*RR2*V2/(E1*RR1*V1+E2*RR2*V2);
-
- Real Zeta = 0;
- if ((rpm1->Zeta!=0)&&(rpm2->Zeta!=0)) { Zeta=(rpm1->Zeta + rpm2->Zeta) / 2.0; }
- if ((rpm1->Zeta!=0)&&(rpm2->Zeta==0)) { Zeta=rpm1->Zeta; }
- if ((rpm1->Zeta==0)&&(rpm2->Zeta!=0)) { Zeta=rpm2->Zeta; }
- if ((rpm1->Zeta==0)&&(rpm2->Zeta==0)) { Zeta=0; }
-
- Real m = 0; //Average mass
- if ((b1->state->mass)&&(b2->state->mass)) {
- m = 0.5*(b1->state->mass+b2->state->mass);
- } else if ((b1->state->mass)&&(!b2->state->mass)) {
- m = b1->state->mass;
- } else if ((!b1->state->mass)&&(b2->state->mass)) {
- m = b2->state->mass;
- }
-
-
- contPhys->cn = 2 * sqrt(contPhys->Kn*m)*Zeta; //Damping coefficient
- //contPhys->cs = contPhys->cn/100.0;
- contPhys->cs = 0.0;
-
-
- if ((rpm1->exampleNumber==rpm2->exampleNumber)&&(initCohesive)&&(scene->iter<=cohesiveThresholdIter)) {
- contPhys->isCohesive=true;
- }
-
- interaction->phys=contPhys;
-}
-
-RpmPhys::~RpmPhys(){};
=== removed file 'pkg/dem/RockPM.hpp'
--- pkg/dem/RockPM.hpp 2012-09-14 20:18:27 +0000
+++ pkg/dem/RockPM.hpp 1970-01-01 00:00:00 +0000
@@ -1,96 +0,0 @@
-/**
-
-=== HIGH LEVEL OVERVIEW OF RockPM ===
-
-Rock Particle Model (RockPM) is a set of classes for modelling
-mechanical behavior of mining rocks.
-*/
-
-#pragma once
-
-#include<yade/pkg/common/Dispatching.hpp>
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/common/PeriodicEngines.hpp>
-#include<yade/pkg/common/NormShearPhys.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-
-class RpmState: public State {
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(RpmState,State,"State information about Rpm body.",
- ((int,specimenNumber,0,,"The variable is used for particle size distribution analyze. Indicates, to which part of specimen belongs para of particles."))
- ((Real,specimenMass,0,,"Indicates the mass of the whole stone, which owns the particle."))
- ((Real,specimenVol,0,,"Indicates the mass of the whole stone, which owns the particle."))
- ((Real,specimenMaxDiam,0,,"Indicates the maximal diametr of the specimen.")),
- /*ctor*/ createIndex();
- );
- REGISTER_CLASS_INDEX(RpmState,State);
-};
-REGISTER_SERIALIZABLE(RpmState);
-
-/** This class holds information associated with each body */
-class RpmMat: public FrictMat {
- public:
- virtual shared_ptr<State> newAssocState() const { return shared_ptr<State>(new RpmState); }
- virtual bool stateTypeOk(State* s) const { return (bool)dynamic_cast<RpmState*>(s); }
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(RpmMat,FrictMat,"Rock material, for use with other Rpm classes.",
- ((int,exampleNumber,0,,"Number of the specimen. This value is equal for all particles of one specimen. [-]"))
- ((bool,initCohesive,false,,"The flag shows, whether particles of this material can be cohesive. [-]"))
- ((Real,stressCompressMax,0,,"Maximal strength for compression. The main destruction parameter. [Pa] //(Needs to be reworked)"))
- ((Real,stressStretchingMax,0,,"Maximal strength for stretching. [Pa]"))
- ((Real,stressShearMax,0,,"Maximal strength for shearing. [Pa]"))
- ((Real,G_over_E,1,,"Ratio of normal/shear stiffness at interaction level. [-]"))
- ((Real,Zeta,0,,"Damping Ratio, http://en.wikipedia.org/wiki/Damping_ratio [-]")),
- createIndex();
- );
-
- REGISTER_CLASS_INDEX(RpmMat,FrictMat);
-};
-REGISTER_SERIALIZABLE(RpmMat);
-
-
-class Ip2_RpmMat_RpmMat_RpmPhys: public IPhysFunctor{
- public:
- virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction);
- FUNCTOR2D(RpmMat,RpmMat);
- DECLARE_LOGGER;
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_RpmMat_RpmMat_RpmPhys,IPhysFunctor,"Convert 2 RpmMat instances to RpmPhys with corresponding parameters.",);
-};
-REGISTER_SERIALIZABLE(Ip2_RpmMat_RpmMat_RpmPhys);
-
-
-class RpmPhys: public NormShearPhys {
- private:
- public:
- Real Fn, epsN, sigmaN, cn, cs;
- Vector3r epsT, Fs, sigmaT;
- bool updatedCnFlag;
-
- virtual ~RpmPhys();
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(RpmPhys,NormShearPhys,"Representation of a single interaction of the Rpm type: storage for relevant parameters.\n\n Evolution of the contact is governed by Law2_Dem3DofGeom_RpmPhys_Rpm, that includes damage effects and chages of parameters inside RpmPhys",
- ((Real,E,NaN,,"normal modulus (stiffness / crossSection) [Pa]"))
- ((Real,crossSection,0,,"equivalent cross-section associated with this contact [m²]"))
- ((Real,G,NaN,,"shear modulus [Pa]"))
- ((Real,tanFrictionAngle,NaN,,"tangens of internal friction angle [-]"))
- ((bool,isCohesive,false,,"if not cohesive, interaction is deleted when distance is greater than lengthMaxTension or less than lengthMaxCompression."))
- ((Real,epsMaxCompression,0,,"Maximal relative penetration of particles during compression. If it is more, the interaction is deleted [m]"))
- ((Real,epsMaxTension,0,,"Maximal relative distance between particles during tension. If it is more, the interaction is deleted [m]"))
- ((Real,epsMaxShear,0,,"Maximal relative distance between particles in tangential direction. If it is more, the interaction is deleted [m]"))
- ((Real,Kn,0,,"Stiffness in normal direction [N/m]"))
- ((Real,Ks,0,,"Stiffness in shear direction [N/m]")),
- /*ctor*/createIndex(); Fn=0; epsN=0; sigmaN=0; sigmaT = Vector3r::Zero(); epsT=Vector3r::Zero(); Fs=Vector3r::Zero(); cn=0; cs=0;
- );
- REGISTER_CLASS_INDEX(RpmPhys,NormShearPhys);
-};
-REGISTER_SERIALIZABLE(RpmPhys);
-
-class Law2_Dem3DofGeom_RockPMPhys_Rpm: public LawFunctor{
- public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
- FUNCTOR2D(Dem3DofGeom,RpmPhys);
-
- YADE_CLASS_BASE_DOC(Law2_Dem3DofGeom_RockPMPhys_Rpm,LawFunctor,"Constitutive law for the Rpm model");
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_RockPMPhys_Rpm);
-
-
=== modified file 'pkg/dem/TriaxialStressController.cpp'
--- pkg/dem/TriaxialStressController.cpp 2013-04-09 11:46:11 +0000
+++ pkg/dem/TriaxialStressController.cpp 2013-07-16 09:57:52 +0000
@@ -23,6 +23,14 @@
Vector3r TriaxialStressController::getStress(int boundId) {assert (boundId>=0 && boundId<=5); return stress[boundId];}
+Vector3r TriaxialStressController::getStrainRate() {
+ return Vector3r (
+ (Body::byId(wall_right_id,scene)->state->vel[0]-Body::byId(wall_left_id,scene)->state->vel[0])/width,
+ (Body::byId(wall_top_id,scene)->state->vel[1]-Body::byId(wall_bottom_id,scene)->state->vel[1])/height,
+ (Body::byId(wall_front_id,scene)->state->vel[2]-Body::byId(wall_back_id,scene)->state->vel[2])/depth
+ );
+}
+
void TriaxialStressController::updateStiffness ()
{
for (int i=0; i<6; ++i) stiffness[i] = 0;
=== modified file 'pkg/dem/TriaxialStressController.hpp'
--- pkg/dem/TriaxialStressController.hpp 2013-02-07 18:12:42 +0000
+++ pkg/dem/TriaxialStressController.hpp 2013-07-16 09:57:52 +0000
@@ -70,8 +70,9 @@
void computeStressStrain();
//! Compute the mean/max unbalanced force in the assembly (normalized by mean contact force)
Real ComputeUnbalancedForce(bool maxUnbalanced=false);
- ///! Getter for stress in python
+ ///! Getter for stress and rates in python
Vector3r getStress(int boundId);
+ Vector3r getStrainRate();
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TriaxialStressController,BoundaryController,
"An engine maintaining constant stresses or constant strain rates on some boundaries of a parallepipedic packing. The stress/strain control is defined for each axis using :yref:`TriaxialStressController::stressMask` (a bitMask) and target values are defined by goal1,goal2, and goal3. sigmaIso has to be defined during growing phases."
@@ -138,6 +139,7 @@
porosity=1;
,
.def_readonly("strain",&TriaxialStressController::strain,"Current strain in a vector (exx,eyy,ezz). The values reflect true (logarithmic) strain.")
+ .def_readonly("strainRate",&TriaxialStressController::getStrainRate,"Current strain rate in a vector d/dt(exx,eyy,ezz).")
.def_readonly("porosity",&TriaxialStressController::porosity,"Porosity of the packing.")
.def_readonly("boxVolume",&TriaxialStressController::boxVolume,"Total packing volume.")
.def_readonly("spheresVolume",&TriaxialStressController::spheresVolume,"Total volume pf spheres.")
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2013-03-07 18:28:01 +0000
+++ pkg/dem/VTKRecorder.cpp 2013-07-23 10:53:20 +0000
@@ -26,7 +26,6 @@
#include<yade/pkg/common/Facet.hpp>
#include<yade/pkg/common/Box.hpp>
#include<yade/pkg/dem/ConcretePM.hpp>
-#include<yade/pkg/dem/RockPM.hpp>
#include<yade/pkg/dem/WirePM.hpp>
#include<yade/pkg/dem/Shop.hpp>
@@ -58,7 +57,6 @@
else if(rec=="mass") recActive[REC_MASS]=true;
else if((rec=="colors") || (rec=="color"))recActive[REC_COLORS]=true;
else if(rec=="cpm") recActive[REC_CPM]=true;
- else if(rec=="rpm") recActive[REC_RPM]=true;
else if(rec=="wpm") recActive[REC_WPM]=true;
else if(rec=="intr") recActive[REC_INTR]=true;
else if((rec=="ids") || (rec=="id")) recActive[REC_ID]=true;
@@ -66,7 +64,7 @@
else if((rec=="clumpids") || (rec=="clumpId")) recActive[REC_CLUMPID]=true;
else if(rec=="materialId") recActive[REC_MATERIALID]=true;
else if(rec=="stress") recActive[REC_STRESS]=true;
- else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, rpm, wpm, intr, id, clumpId, materialId). Ignored.");
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId). Ignored.");
}
// cpm needs interactions
if(recActive[REC_CPM]) recActive[REC_INTR]=true;
@@ -200,17 +198,6 @@
cpmStress->SetNumberOfComponents(9);
cpmStress->SetName("cpmStress");
- // extras for RPM
- vtkSmartPointer<vtkFloatArray> rpmSpecNum = vtkSmartPointer<vtkFloatArray>::New();
- rpmSpecNum->SetNumberOfComponents(1);
- rpmSpecNum->SetName("rpmSpecNum");
- vtkSmartPointer<vtkFloatArray> rpmSpecMass = vtkSmartPointer<vtkFloatArray>::New();
- rpmSpecMass->SetNumberOfComponents(1);
- rpmSpecMass->SetName("rpmSpecMass");
- vtkSmartPointer<vtkFloatArray> rpmSpecDiam = vtkSmartPointer<vtkFloatArray>::New();
- rpmSpecDiam->SetNumberOfComponents(1);
- rpmSpecDiam->SetName("rpmSpecDiam");
-
// extras for WireMatPM
vtkSmartPointer<vtkFloatArray> wpmNormalForce = vtkSmartPointer<vtkFloatArray>::New();
wpmNormalForce->SetNumberOfComponents(1);
@@ -356,11 +343,6 @@
float s[9]={ (float) ss(0,0), (float) ss(0,1), (float) ss(0,2), (float) ss(1,0), (float) ss(1,1), (float) ss(1,2), (float) ss(2,0), (float) ss(2,1), (float) ss(2,2)};
cpmStress->InsertNextTupleValue(s);
}
- if (recActive[REC_RPM]){
- rpmSpecNum->InsertNextValue(YADE_PTR_CAST<RpmState>(b->state)->specimenNumber);
- rpmSpecMass->InsertNextValue(YADE_PTR_CAST<RpmState>(b->state)->specimenMass);
- rpmSpecDiam->InsertNextValue(YADE_PTR_CAST<RpmState>(b->state)->specimenMaxDiam);
- }
if (recActive[REC_MATERIALID]) spheresMaterialId->InsertNextValue(b->material->id);
continue;
@@ -480,11 +462,6 @@
spheresUg->GetPointData()->AddArray(cpmDamage);
spheresUg->GetPointData()->AddArray(cpmStress);
}
- if (recActive[REC_RPM]){
- spheresUg->GetPointData()->AddArray(rpmSpecNum);
- spheresUg->GetPointData()->AddArray(rpmSpecMass);
- spheresUg->GetPointData()->AddArray(rpmSpecDiam);
- }
if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
=== modified file 'py/export.py'
--- py/export.py 2013-04-08 15:51:39 +0000
+++ py/export.py 2013-07-10 09:16:15 +0000
@@ -163,6 +163,8 @@
outFile.close()
self.snapCount+=1
+
+#text===============================================================
def text(filename,mask=-1):
"""Save sphere coordinates into a text file; the format of the line is: x y z r. Non-spherical bodies are silently skipped. Example added to examples/regular-sphere-pack/regular-sphere-pack.py
@@ -192,12 +194,13 @@
def __init__(self,baseName,startSnap=0):
self.spheresSnapCount = startSnap
self.facetsSnapCount = startSnap
+ self.intrsSnapCount = startSnap
self.baseName = baseName
def exportSpheres(self,ids='all',what=[],comment="comment",numLabel=None):
"""exports spheres (positions and radius) and defined properties.
- :param ids: if dd "all", then export all spheres, otherwise only spheres from integer list
+ :param ids: if "all", then export all spheres, otherwise only spheres from integer list
:type ids: [int] | "all"
:param what: what other than then position and radius export. parameter is list of couple (name,command). Name is string under which it is save to vtk, command is string to evaluate. Node that the bodies are labeled as b in this function. Scalar, vector and tensor variables are supported. For example, to export velocity (with name particleVelocity) and the distance form point (0,0,0) (named as dist) you should write: ... what=[('particleVelocity','b.state.vel'),('dist','b.state.pos.norm()', ...
:type what: [tuple(2)]
@@ -378,6 +381,66 @@
outFile.close()
self.facetsSnapCount += 1
+ def exportInteractions(self,ids='all',what=[],comment="comment",numLabel=None):
+ """exports interactions and defined properties.
+
+ :param ids: if "all", then export all interactions, otherwise only interactions from (integer,integer) list
+ :type ids: [(int,int)] | "all"
+ :param what: what to export. parameter is list of couple (name,command). Name is string under which it is save to vtk, command is string to evaluate. Node that the interactions are labeled as i in this function. Scalar, vector and tensor variables are supported. For example, to export stiffness difference from certain value (1e9) (named as dStiff) you should write: ... what=[('dStiff','i.phys.kn-1e9'), ...
+ :type what: [tuple(2)]
+ :param string comment: comment to add to vtk file
+ :param int numLabel: number of file (e.g. time step), if unspecified, the last used value + 1 will be used
+ """
+ allIds = False
+ if ids=='all':
+ intrs=O.interactions
+ allIds = True
+ else:
+ intrs = []
+ for j,k in ids:
+ i = O.interactions[j,k]
+ if not i: continue
+ intrs.append(i)
+ n = len([i for i in intrs if i.isReal])
+ ids = [None for j in xrange(n)]
+ outFile = open(self.baseName+'-intrs-%04d'%self.intrsSnapCount+'.vtk', 'w')
+ outFile.write("# vtk DataFile Version 3.0.\n%s\nASCII\n\nDATASET POLYDATA\nPOINTS %d double\n"%(comment,2*n))
+ for j,i in enumerate(intrs):
+ if not i.isReal: continue
+ pos = O.bodies[i.id1].state.pos
+ outFile.write("%g %g %g\n"%(pos[0],pos[1],pos[2]))
+ pos = O.bodies[i.id2].state.pos + i.cellDist
+ outFile.write("%g %g %g\n"%(pos[0],pos[1],pos[2]))
+ ids[j] = (2*j-1,2*j)
+ outFile.write("LINES %d %d\n"%(n,3*n))
+ for j,i in enumerate(intrs):
+ outFile.write("2 %d %d\n"%(ids[j][0]+1,ids[j][1]+1))
+ outFile.write("\nCELL_DATA %d\n"%(n))
+ for name,command in what:
+ test = eval(command)
+ if isinstance(test,Matrix3):
+ print 'WARNING: export.VTKExporter.exportInteractions: wrong `what` parameter, Matrix3 output not (yet?) supported'
+ #outFile.write("\nTENSORS %s double\n"%(name))
+ #for i in intrs:
+ # t = eval(command)
+ # outFile.write("%g %g %g\n%g %g %g\n%g %g %g\n\n"%(t[0,0],t[0,1],t[0,2],t[1,0],t[1,1],t[1,2],t[2,0],t[2,1],t[2,2]))
+ elif isinstance(test,Vector3):
+ print 'WARNING: export.VTKExporter.exportInteractions: wrong `what` parameter, Vector3 output not (yet?) supported'
+ #outFile.write("\nVECTORS %s double\n"%(name))
+ #for i in intrs:
+ # v = eval(command)
+ # outFile.write("%g %g %g\n"%(v[0],v[1],v[2]))
+ elif isinstance(test,(int,float)):
+ outFile.write("\nSCALARS %s double 1\nLOOKUP_TABLE default\n"%(name))
+ for i in intrs:
+ outFile.write("%g\n"%(eval(command)))
+ else:
+ print 'WARNING: export.VTKExporter.exportInteractions: wrong `what` parameter, vtk output might be corrupted'
+ outFile.close()
+ self.intrsSnapCount += 1
+
+
+
#gmshGeoExport===============================================================
def gmshGeo(filename, comment='',mask=-1,accuracy=-1):
"""Save spheres in geo-file for the following using in GMSH (http://www.geuz.org/gmsh/doc/texinfo/) program. The spheres can be there meshed.
=== modified file 'py/mathWrap/miniEigen.cpp'
--- py/mathWrap/miniEigen.cpp 2013-03-07 18:28:01 +0000
+++ py/mathWrap/miniEigen.cpp 2013-07-09 12:42:00 +0000
@@ -125,8 +125,10 @@
#include<Eigen/SVD>
static bp::tuple Matrix3r_polarDecomposition(const Matrix3r& self){ Matrix3r unitary,positive; Matrix_computeUnitaryPositive(self,&unitary,&positive); return bp::make_tuple(unitary,positive); }
static bp::tuple Matrix3r_eigenDecomposition(const Matrix3r& self){ Matrix3r rot,diag; matrixEigenDecomposition(self,rot,diag); return bp::make_tuple(rot,Vector3r(diag.diagonal())); }
+static bp::tuple Matrix3r_svd(const Matrix3r& self){ Matrix3r u,s,v; Matrix_SVD(self,&u,&s,&v); return bp::make_tuple(u,s,v); }
static bp::tuple Matrix6r_polarDecomposition(const Matrix6r& self){ Matrix6r unitary,positive; Matrix_computeUnitaryPositive(self,&unitary,&positive); return bp::make_tuple(unitary,positive); }
static bp::tuple Matrix6r_eigenDecomposition(const Matrix6r& self){ Matrix6r rot,diag; matrixEigenDecomposition(self,rot,diag); return bp::make_tuple(rot,Vector6r(diag.diagonal())); }
+static bp::tuple Matrix6r_svd(const Matrix6r& self){ Matrix6r u,s,v; Matrix_SVD(self,&u,&s,&v); return bp::make_tuple(u,s,v); }
#undef IDX_CHECK
@@ -237,6 +239,7 @@
.def("transpose",&Matrix3r_transpose)
.def("polarDecomposition",&Matrix3r_polarDecomposition)
.def("eigenDecomposition",&Matrix3r_eigenDecomposition)
+ .def("svd",&Matrix3r_svd)
.def("diagonal",&Matrix3r_diagonal)
.def("row",&Matrix3r_row)
.def("col",&Matrix3r_col)
@@ -275,6 +278,7 @@
.def("col",&Matrix6r_col)
.def("polarDecomposition",&Matrix6r_polarDecomposition)
.def("eigenDecomposition",&Matrix6r_eigenDecomposition)
+ .def("svd",&Matrix6r_svd)
//
.def("__neg__",&Matrix6r__neg__)
.def("__add__",&Matrix6r__add__Matrix6r).def("__iadd__",&Matrix6r__iadd__Matrix6r)
=== modified file 'py/tests/clump.py'
--- py/tests/clump.py 2013-03-08 21:26:47 +0000
+++ py/tests/clump.py 2013-07-08 08:33:16 +0000
@@ -43,7 +43,9 @@
# centroid
S=b1.state.mass*b1.state.pos+b2.state.mass*b2.state.pos
c=S/bC.state.mass
- self.assertEqual(bC.state.pos,c);
+ self.assertAlmostEqual(bC.state.pos[0],c[0]);
+ self.assertAlmostEqual(bC.state.pos[1],c[1]);
+ self.assertAlmostEqual(bC.state.pos[2],c[2]);
# inertia
i1,i2=(8./15)*pi*b1.material.density*b1.shape.radius**5, (8./15)*pi*b2.material.density*b2.shape.radius**5 # inertia of spheres
iMax=i1+i2+b1.state.mass*(b1.state.pos-c).norm()**2+b2.state.mass*(b2.state.pos-c).norm()**2 # minimum principal inertia
=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp 2013-06-21 06:34:26 +0000
+++ py/wrapper/yadeWrapper.cpp 2013-07-23 10:01:46 +0000
@@ -647,6 +647,19 @@
shared_ptr<Scene> scene_get(){ return OMEGA.getScene(); }
int addScene(){return OMEGA.addScene();}
void switchToScene(int i){OMEGA.switchToScene(i);}
+ string sceneToString(){
+ ostringstream oss;
+ yade::ObjectIO::save<typeof(OMEGA.getScene()),boost::archive::binary_oarchive>(oss,"scene",OMEGA.getScene());
+ oss.flush();
+ return oss.str();
+ }
+ void stringToScene(const string &sstring, string mark=""){
+ Py_BEGIN_ALLOW_THREADS; OMEGA.stop(); Py_END_ALLOW_THREADS;
+ assertScene();
+ OMEGA.memSavedSimulations[":memory:"+mark]=sstring;
+ OMEGA.sceneFile=":memory:"+mark;
+ load(OMEGA.sceneFile,true);
+ }
void save(std::string fileName,bool quiet=false){
assertScene();
@@ -802,6 +815,8 @@
.def("addScene",&pyOmega::addScene,"Add new scene to Omega, returns its number")
.def("switchToScene",&pyOmega::switchToScene,"Switch to defined scene. Default scene has number 0, other scenes have to be created by addScene method.")
.def("switchScene",&pyOmega::switchScene,"Switch to alternative simulation (while keeping the old one). Calling the function again switches back to the first one. Note that most variables from the first simulation will still refer to the first simulation even after the switch\n(e.g. b=O.bodies[4]; O.switchScene(); [b still refers to the body in the first simulation here])")
+ .def("sceneToString",&pyOmega::sceneToString,"Return the entire scene as a string. Equivalent to using O.save(...) except that the scene goes to a string instead of a file. (see also stringToScene())")
+ .def("stringToScene",&pyOmega::stringToScene,(python::arg("mark")=""),"Load simulation from a string passed as argument (see also sceneToString).")
.def("labeledEngine",&pyOmega::labeled_engine_get,"Return instance of engine/functor with the given label. This function shouldn't be called by the user directly; every ehange in O.engines will assign respective global python variables according to labels.\n\nFor example:\n\n\t *O.engines=[InsertionSortCollider(label='collider')]*\n\n\t *collider.nBins=5 # collider has become a variable after assignment to O.engines automatically*")
.def("resetTime",&pyOmega::resetTime,"Reset simulation time: step number, virtual and real time. (Doesn't touch anything else, including timings).")
.def("plugins",&pyOmega::plugins_get,"Return list of all plugins registered in the class factory.")
