yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #09942
Re: [Branch ~yade-pkg/yade/git-trunk] Rev 3687: Merge branch 'master' of github.com:yade/trunk
Hello,
is this ok, or did I do something wrong? the same question for the next
noreply message..
Thanks
Jan
2013/9/8 <noreply@xxxxxxxxxxxxx>
> Merge authors:
> Anton Gladky (gladky-anton)
> Bruno Chareyre (bruno-chareyre)
> Klaus Thoeni (klaus.thoeni)
> ------------------------------------------------------------
> revno: 3687 [merge]
> committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
> timestamp: Sun 2013-09-08 13:19:06 +0200
> message:
> Merge branch 'master' of github.com:yade/trunk
> removed:
> examples/test/CundallStrackTest.py
> examples/test/Dem3DofGeom.py
> pkg/dem/Dem3DofGeom_FacetSphere.cpp
> pkg/dem/Dem3DofGeom_FacetSphere.hpp
> pkg/dem/Dem3DofGeom_SphereSphere.cpp
> pkg/dem/Dem3DofGeom_SphereSphere.hpp
> pkg/dem/Dem3DofGeom_WallSphere.cpp
> pkg/dem/Dem3DofGeom_WallSphere.hpp
> py/_eudoxos.cpp
> py/eudoxos.py
> added:
> scripts/checks-and-tests/checks/DEM-PFV-check.py
> scripts/checks-and-tests/checks/data/100spheres
> modified:
> CMakeLists.txt
> cMake/FindMetis.cmake
> core/Clump.cpp
> core/ForceContainer.hpp
> core/main/main.py.in
> doc/references.bib
> doc/sphinx/formulation.rst
> doc/sphinx/introduction.rst
> doc/sphinx/prog.rst
> doc/sphinx/references.bib
> doc/sphinx/templates/index.html
> doc/sphinx/user.rst
> doc/yade-articles.bib
> doc/yade-conferences.bib
> doc/yade-theses.bib
> pkg/common/Cylinder.hpp
> pkg/dem/ConcretePM.cpp
> pkg/dem/DemXDofGeom.cpp
> pkg/dem/DemXDofGeom.hpp
> pkg/dem/DomainLimiter.cpp
> pkg/dem/FlowEngine.cpp
> pkg/dem/FlowEngine.hpp
> pkg/dem/HertzMindlin.cpp
> pkg/dem/HertzMindlin.hpp
> pkg/dem/Shop.cpp
> pkg/dem/Shop.hpp
> pkg/dem/VTKRecorder.hpp
> pkg/dem/ViscoelasticCapillarPM.cpp
> pkg/dem/ViscoelasticPM.cpp
> pkg/dem/ViscoelasticPM.hpp
> py/CMakeLists.txt
> py/_extraDocs.py
> py/_utils.cpp
> py/system.py
> py/tests/__init__.py
> scripts/checks-and-tests/checks/README
> scripts/checks-and-tests/checks/checkGravity.py
> scripts/checks-and-tests/checks/checkList.py
> scripts/checks-and-tests/checks/checkTestDummy.py
> scripts/checks-and-tests/checks/checkTestTriax.py
> scripts/checks-and-tests/checks/checkWirePM.py
>
>
> --
> 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-08-15 12:52:44 +0000
> +++ CMakeLists.txt 2013-09-05 07:52:37 +0000
> @@ -380,7 +380,6 @@
> SET (pyExecutable ${PYTHON_EXECUTABLE})
> SET (profile "default")
> SET (sourceRoot "${CMAKE_CURRENT_SOURCE_DIR}")
> -SET (libstdcxx "/usr/lib/libstdc++.so.6") #This is a hack, which is
> not needed for modern systems. Should be deleted.
> #====================================
>
> CONFIGURE_FILE(core/main/yade-batch.in"${CMAKE_BINARY_DIR}/bins/yade${SUFFIX}-batch")
> @@ -409,6 +408,7 @@
> MESSAGE(STATUS "Disabled features:${DISABLED_FEATS}")
> IF (DEBUG)
> MESSAGE(STATUS "Debug build")
> + SET (debugbuild " (debug build)")
> ELSE (DEBUG)
> MESSAGE(STATUS "Optimized build")
> ENDIF (DEBUG)
>
> === modified file 'cMake/FindMetis.cmake'
> --- cMake/FindMetis.cmake 2013-06-26 18:15:55 +0000
> +++ cMake/FindMetis.cmake 2013-08-30 13:49:57 +0000
> @@ -5,8 +5,8 @@
> # METIS_LIBRARY, libraries to link against to use GL2PS.
> # METIS_FOUND, If false, do not try to use GL2PS.
>
> -FIND_PATH(METIS_INCLUDE_DIR metis.h)
> -FIND_LIBRARY(METIS_LIBRARY NAMES metis )
> +FIND_PATH(METIS_INCLUDE_DIR metis.h PATHS /usr/include/metis)
> +FIND_LIBRARY(METIS_LIBRARY NAMES metis PATHS /usr/lib)
>
> # handle the QUIETLY and REQUIRED arguments and set LOKI_FOUND to TRUE if
> # all listed variables are TRUE
>
> === modified file 'core/Clump.cpp'
> --- core/Clump.cpp 2013-07-05 07:33:47 +0000
> +++ core/Clump.cpp 2013-08-30 06:19:46 +0000
> @@ -119,7 +119,7 @@
> const Sphere* sphere1 = YADE_CAST<Sphere*>
> (subBody1->shape.get());
> const Sphere* sphere2 = YADE_CAST<Sphere*>
> (subBody2->shape.get());
> Real un =
> (sphere1->radius+sphere2->radius) - dist.norm();
> - if (un > 0.) {intersecting = true; break;}
> + if (un >
> -0.001*min(sphere1->radius,sphere2->radius)) {intersecting = true; break;}
> }
> }
> }
>
> === modified file 'core/ForceContainer.hpp'
> --- core/ForceContainer.hpp 2013-08-22 18:35:55 +0000
> +++ core/ForceContainer.hpp 2013-08-26 21:41:10 +0000
> @@ -72,7 +72,7 @@
> // dummy function to avoid template resolution failure
> friend class boost::serialization::access; template<class
> ArchiveT> void serialize(ArchiveT & ar, unsigned int version){}
> public:
> - ForceContainer(): size(0),
> permSize(0),syncedSizes(true),synced(true),moveRotUsed(false),_zero(Vector3r::Zero()),syncCount(0),lastReset(0){
> + ForceContainer(): size(0),
> permSize(0),syncedSizes(true),synced(true),moveRotUsed(false),permForceUsed(false),_zero(Vector3r::Zero()),syncCount(0),lastReset(0){
> nThreads=omp_get_max_threads();
> for(int i=0; i<nThreads; i++){
> _forceData.push_back(vvector());
> _torqueData.push_back(vvector());
> @@ -220,7 +220,7 @@
> // dummy function to avoid template resolution failure
> friend class boost::serialization::access; template<class
> ArchiveT> void serialize(ArchiveT & ar, unsigned int version){}
> public:
> - ForceContainer(): size(0), permSize(0),
> moveRotUsed(false), syncCount(0), lastReset(0){}
> + ForceContainer(): size(0), permSize(0),
> moveRotUsed(false), permForceUsed(false), syncCount(0), lastReset(0){}
> const Vector3r& getForce(Body::id_t id){ensureSize(id);
> return _force[id];}
> void addForce(Body::id_t id,const Vector3r&
> f){ensureSize(id); _force[id]+=f;}
> const Vector3r& getTorque(Body::id_t id){ensureSize(id);
> return _torque[id];}
>
> === modified file 'core/main/main.py.in'
> --- core/main/main.py.in 2013-05-07 09:21:35 +0000
> +++ core/main/main.py.in 2013-09-05 07:52:37 +0000
> @@ -11,7 +11,7 @@
> # get yade path (allow YADE_PREFIX to override)
> prefix,suffix='${runtimePREFIX}' if not os.environ.has_key('YADE_PREFIX')
> else os.environ['YADE_PREFIX'],'${SUFFIX}'
> # duplicate some items from yade.config here, so that we can increase
> verbosity when the c++ part is booting
> -features,version='${features}'.split(' '),'${realVersion}'
> +features,version,debugbuild='${features}'.split(' '),'${realVersion}','
> ${debugbuild}'
>
> libPATH='${LIBRARY_OUTPUT_PATH}'
> if (libPATH[1:] == '{LIBRARY_OUTPUT_PATH}'): libPATH='lib'
> @@ -52,11 +52,17 @@
> par.add_argument('--performance',help='Starts a test to measure the
> productivity',dest='performance',action='store_true')
>
> par.add_argument('script',nargs='?',default='',type=str,help=argparse.SUPPRESS)
> par.add_argument('args',nargs=argparse.REMAINDER,help=argparse.SUPPRESS)
> # see argparse doc, par.disable_interspersed_args() from optargs module
> +par.add_argument('-l',help='import libraries at startup before importing
> yade libs. May be used when the ordering of imports matter (see e.g.
> https://bugs.launchpad.net/yade/+bug/1183402/comments/3). The option can
> be use multiple times, as in "yade -llib1
> -llib2"',default=None,action='append',dest='impLibraries',type=str)
> opts=par.parse_args()
> args = opts.args
>
> +if opts.impLibraries:
> + sys.path.append('.')
> + for lib in opts.impLibraries:
> + __import__(lib)
> +
> if opts.version:
> - print 'Yade version: %s'%version
> + print 'Yade version: %s%s'%(version,debugbuild)
> sys.exit(0)
>
> if opts.script:
> @@ -109,7 +115,7 @@
> else: os.environ['OMP_NUM_THREADS']='1'
>
> if __name__ == "__main__": # do not print this while importing yade in
> other python application
> - sys.stderr.write('Welcome to Yade '+version+'\n')
> + sys.stderr.write('Welcome to Yade '+version+debugbuild+'\n')
>
> # initialization and c++ plugins import
> import yade
>
> === modified file 'doc/references.bib'
> --- doc/references.bib 2013-08-14 08:01:55 +0000
> +++ doc/references.bib 2013-08-28 10:26:24 +0000
> @@ -580,18 +580,6 @@
> url = "http://www.refdoc.fr/Detailnotice?idarticle=7435486";
> }
>
> -@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{Luding2008,
> year={2008},
> issn={1434-5021},
>
> === modified file 'doc/sphinx/formulation.rst'
> --- doc/sphinx/formulation.rst 2013-05-14 21:24:11 +0000
> +++ doc/sphinx/formulation.rst 2013-08-29 10:30:31 +0000
> @@ -314,8 +314,6 @@
> -------------
> In order to keep $\vec{u}_T$ consistent (e.g. that $\vec{u}_T$ must be
> constant if two spheres retain mutually constant configuration but move
> arbitrarily in space), then either $\vec{u}_T$ must track spheres' spatial
> motion or must (somehow) rely on sphere-local data exclusively.
>
> -These two possibilities lead to two algorithms of computing shear
> strains. They should give the same results (disregarding numerical
> imprecision), but there is a trade-off between computational cost of the
> incremental method and robustness of the total one.
> -
> Geometrical meaning of shear strain is shown in `fig-shear-2d`_.
>
> .. _fig-shear-2d:
> @@ -323,9 +321,7 @@
>
> Evolution of shear displacement $\vec{u}_T$ due to mutual motion
> of spheres, both linear and rotational. Left configuration is the initial
> contact, right configuration is after displacement and rotation of one
> particle.
>
> -Incremental algorithm
> -^^^^^^^^^^^^^^^^^^^^^
> -The incremental algorithm is widely used in DEM codes and is described
> frequently ([Luding2008]_, [Alonso2004]_). Yade implements this algorithm
> in the :yref:`ScGeom` class. At each step, shear displacement $\uT$ is
> updated; the update increment can be decomposed in 2 parts: motion of the
> interaction (i.e. $\vec{C}$ and $\vec{n}$) in global space and mutual
> motion of spheres.
> +The classical incremental algorithm is widely used in DEM codes and is
> described frequently ([Luding2008]_, [Alonso2004]_). Yade implements this
> algorithm in the :yref:`ScGeom` class. At each step, shear displacement
> $\uT$ is updated; the update increment can be decomposed in 2 parts: motion
> of the interaction (i.e. $\vec{C}$ and $\vec{n}$) in global space and
> mutual motion of spheres.
>
> #. Contact moves dues to changes of the spheres' positions $\vec{C}_1$
> and $\vec{C}_2$, which updates current $\currC$ and $\currn$ as per
> :eq:`eq-contact-point` and :eq:`eq-contact-normal`. $\prevuT$ is
> perpendicular to the contact plane at the previous step $\prevn$ and must
> be updated so that $\prevuT+(\Delta\uT)=\curruT\perp\currn$; this is done
> by perpendicular projection to the plane first (which might decrease
> $|\uT|$) and adding what corresponds to spatial rotation of the interaction
> instead:
>
> @@ -353,74 +349,6 @@
>
> .. math:: \curruT=\prevuT+(\Delta\uT)_1 + (\Delta\uT)_2 + (\Delta\uT)_3.
>
> -.. _sect-formulation-total-shear:
> -
> -Total algorithm
> -^^^^^^^^^^^^^^^
> -The following algorithm, aiming at stabilization of response even with
> large rotation speeds or $\Delta t$ approaching stability limit, was
> designed in [Smilauer2010b]_. (A similar algorithm based on total
> formulation, which covers additionally bending and torsion, was proposed in
> [Wang2009]_.) It is based on tracking original contact points (with zero
> shear) in the particle-local frame.
> -
> -In this section, variable symbols implicitly denote their current values
> unless explicitly stated otherwise.
> -
> -Shear strain may have two sources: mutual rotation of spheres or
> transversal displacement of one sphere with respect to the other. Shear
> strain does not change if both spheres move or rotate but are not in linear
> or angular motion mutually. To accurately and reliably model this
> situation, for every new contact the initial contact point $\bar{\vec{C}}$
> is mapped into local sphere coordinates ($\vec{p}_{01}$, $\vec{p}_{02}$).
> As we want to determine the distance between both points (i.e. how long the
> trajectory in on both spheres' surfaces together), the shortest path from
> current $\vec{C}$ to the initial locally mapped point on the sphere's
> surface is „unrolled“ to the contact plane ($\vec{p}'_{01}$,
> $\vec{p}'_{02}$); then we can measure their linear distance $\uT$ and
> define shear strain $\vec{\eps}_T=\uT/d_0$ (fig. `fig-shear-displacement`_).
> -
> -More formally, taking $\bar{\vec{C}}_i$, $\bar{q}_i$ for the sphere
> initial positions and orientations (as quaterions) in global coordinates,
> the initial sphere-local contact point *orientation* (relative to
> sphere-local axis $\hat{x}$) is remembered:
> -
> -.. math::
> - :nowrap:
> -
> - \begin{align*}
> - \bar{\vec{n}}&=\normalized{{\vec{C}}_1-{\vec{C}}_2}, \\
> - \bar{q}_{01}&=\Align(\hat x,\bar{q}_1^*\bar{\vec{n}}
> \bar{q}_1^{**}), \\
> - \bar{q}_{02}&=\Align(\hat x,\bar{q}_2^* (-\bar{\vec{n}})
> \bar{q}_2^{**}).
> - \end{align*}
> -
> -.. (See \autoref{sect-quaternions} for definition of $\Align$.)
> -
> -
> -After some spheres motion, the original point can be "unrolled" to the
> current contact plane:
> -
> -.. math::
> - :nowrap:
> -
> - \begin{align*}
> - q&=\Align(\vec{n},q_1 \bar{q}_{01} \hat x (q_1
> \bar{q}_{01})^*) \quad\hbox{(auxiliary)} \\
> - \vec{p}'_{01}&=q_{\theta}d_1(q_{\vec{u}} \times \vec{n})
> - \end{align*}
> -
> -where $q_{\vec{u}}$, $q_{\theta}$ are axis and angle components of $q$
> and $p_{01}'$ is the unrolled point. Similarly,
> -
> -.. math::
> - :nowrap:
> -
> - \begin{align*}
> - q&=\Align(\vec{n},q_2 \bar{q}_{02} \hat x (q_2
> \bar{q}_{02})^*) \\
> - \vec{p}'_{02}&=q_{\theta}d_1(q_{\vec{u}} \times
> (-\vec{n})).
> - \end{align*}
> -
> -Shear displacement and strain are then computed easily:
> -
> -.. math::
> - :nowrap:
> -
> - \begin{align*}
> - \uT&=\vec{p}'_{02}-\vec{p}'_{01} \\
> - \vec{\eps}_T&=\frac{\uT}{d_0}
> - \end{align*}
> -
> -When using material law with plasticity in shear, it may be necessary to
> limit maximum shear strain, in which case the mapped points are moved
> closer together to the requested distance (without changing
> $\hat{\vec{u}}_T$). This allows us to remember the previous strain
> direction and also avoids summation of increments of plastic strain at
> every step (`fig-shear-slip`_).
> -
> -.. _fig-shear-displacement:
> -.. figure:: fig/shear-displacement.*
> -
> - Shear displacement computation for two spheres in relative motion.
> -
> -
> -.. _fig-shear-slip:
> -.. figure:: fig/shear-slip.*
> -
> - Shear plastic slip for two spheres.
> -
> -This algorithm is straightforwardly modified to facet-sphere
> interactions. In Yade, it is implemented by :yref:`Dem3DofGeom` and related
> classes.
>
> .. _sect-formulation-stress-cundall:
>
>
> === modified file 'doc/sphinx/introduction.rst'
> --- doc/sphinx/introduction.rst 2013-08-14 08:01:55 +0000
> +++ doc/sphinx/introduction.rst 2013-08-29 10:30:31 +0000
> @@ -343,7 +343,7 @@
> :yref:`IGeom`
> holding geometrical configuration of the two particles in
> collision; it is updated automatically as the particles in question move
> and can be queried for various geometrical characteristics, such as
> penetration distance or shear strain.
>
> - Based on combination of types of :yref:`Shapes<Shape>` of the
> particles, there might be different storage requirements; for that reason,
> a number of derived classes exists, e.g. for representing geometry of
> contact between :yref:`Sphere+Sphere<Dem3DofGeom_SphereSphere>`,
> :yref:`Facet+Sphere<Dem3DofGeom_FacetSphere>` etc.
> + Based on combination of types of :yref:`Shapes<Shape>` of the
> particles, there might be different storage requirements; for that reason,
> a number of derived classes exists, e.g. for representing geometry of
> contact between :yref:`Sphere+Sphere<ScGeom>`,
> :yref:`Cylinder+Sphere<CylScGeom>` etc. Note, however, that it is possible
> to represent many type of contacts with the basic sphere-sphere geometry
> (for instance in :yref:`Ig2_Wall_Sphere_ScGeom`).
> :yref:`IPhys`
> representing non-geometrical features of the interaction; some are
> computed from :yref:`Materials<Material>` of the particles in contact using
> some averaging algorithm (such as contact stiffness from Young's moduli of
> particles), others might be internal variables like damage.
>
>
> === modified file 'doc/sphinx/prog.rst'
> --- doc/sphinx/prog.rst 2013-08-14 08:01:55 +0000
> +++ doc/sphinx/prog.rst 2013-08-29 10:30:31 +0000
> @@ -872,7 +872,7 @@
>
> .. code-block:: c++
>
> - class Ig2_Facet_Sphere_Dem3DofGeom: public IGeomFunctor{
> + class Ig2_Facet_Sphere_ScGeom: public IGeomFunctor{
> public:
>
> // override the IGeomFunctor::go
> @@ -914,7 +914,7 @@
>
> If no functor is able to accept given types (first condition violated) or
> multiple functors have the same distance (in condition 2), an exception is
> thrown.
>
> -This resolution mechanism makes it possible, for instance, to have a
> hierarchy of :yref:`Dem3DofGeom` classes (for different combination of
> shapes: :yref:`Dem3DofGeom_SphereSphere`, :yref:`Dem3DofGeom_FacetSphere`,
> :yref:`Dem3DofGeom_WallSphere`), but only provide a :yref:`LawFunctor`
> accepting ``Dem3DofGeom``, rather than having different laws for each shape
> combination.
> +This resolution mechanism makes it possible, for instance, to have a
> hierarchy of :yref:`ScGeom` classes (for different combination of shapes),
> but only provide a :yref:`LawFunctor` accepting ``ScGeom``, rather than
> having different laws for each shape combination.
>
> .. note:: Performance implications of dispatch resolution are relatively
> low. The dispatcher lookup is only done once, and uses fast lookup matrix
> (1D or 2D); then, the functor found for this type(s) is cached within the
> ``Interaction`` (or ``Body``) instance. Thus, regular functor call costs
> the same as dereferencing pointer and calling virtual method. There is
> `blueprint <
> https://blueprints.launchpad.net/yade/+spec/devirtualize-functor-calls>`__
> to avoid virtual function call as well.
>
>
> === modified file 'doc/sphinx/references.bib'
> --- doc/sphinx/references.bib 2013-06-06 09:27:41 +0000
> +++ doc/sphinx/references.bib 2013-08-28 10:38:02 +0000
> @@ -305,12 +305,5 @@
> doi = "10.1016/S0378-4371(99)00183-1",
> url = "
> http://www.sciencedirect.com/science/article/pii/S0378437199001831";,
> author = "Y.C. Zhou and B.D. Wright and R.Y. Yang and B.H. Xu and
> A.B. Yu",
> - keywords = "Static sandpiles",
> - keywords = "Granular compaction",
> - keywords = "Dynamics and kinematics of a particle and a system of
> particles",
> - keywords = "Granular flow",
> - keywords = "Mixing",
> - keywords = "Segregation and stratification",
> - keywords = "Granular systems "
> }
>
>
> === modified file 'doc/sphinx/templates/index.html'
> --- doc/sphinx/templates/index.html 2013-02-26 19:28:47 +0000
> +++ doc/sphinx/templates/index.html 2013-08-30 19:39:27 +0000
> @@ -8,6 +8,12 @@
> Yade is located at <a href="https://www.yade-dem.org";>www.yade-dem.org</a>,
> which contains <a href="https://www.yade-dem.org/doc/";>this
> documentation</a> and <a href="https://www.yade-dem.org/wiki/";>wiki</a>.
> Development is kindly hosted at <a href="http://www.launchpad.net/yade";>launchpad</a>;
> it is used for <a href="http://code.launchpad.net/yade";>source code</a>,
> <a href="http://bugs.launchpad.net/yade";>bug tracking</a>, <a href="
> http://blueprints.launchpad.net/yade";>planning</a>, <a href="
> https://launchpad.net/~yade-pkg/+archive/snapshots";>package</a> and <a
> href="https://launchpad.net/yade/+download";>source</a> downloads</a> and
> more.
> </p>
> <p>
> + Since March 2012 Yade is using GIT as VCS and placing the source code on
> + <a href="https://github.com/yade/trunk";>github</a>. The Launchpad is
> importing
> + source code several times per day for the further buiding of
> "daily"-packages
> + for Ubuntu.
> + </p>
> + <p>
> Please make sure you read the <a href="{{ pathto("citing")
> }}">"Acknowledging Yade"</a> section if you plan to cite Yade in
> publications.
> </p>
> <p>
>
> === modified file 'doc/sphinx/user.rst'
> --- doc/sphinx/user.rst 2013-08-22 18:35:55 +0000
> +++ doc/sphinx/user.rst 2013-08-29 10:30:31 +0000
> @@ -418,7 +418,7 @@
> #. Approximate collision detection is adjusted so that approximate
> contacts are detected also between particles within the interaction radius.
> This consists in setting value of :yref:`Bo1_Sphere_Aabb.aabbEnlargeFactor`
> to the interaction radius value.
>
> #. The geometry functor (``Ig2``)
> - would normally say that "there is no contact" if given 2 spheres that
> are not in contact. Therefore, the same value as for
> :yref:`Bo1_Sphere_Aabb.aabbEnlargeFactor` must be given to it. (Either
> :yref:`Ig2_Sphere_Sphere_Dem3DofGeom.distFactor` or
> :yref:`Ig2_Sphere_Sphere_ScGeom.interactionDetectionFactor`, depending on
> the functor that is in use.
> + would normally say that "there is no contact" if given 2 spheres that
> are not in contact. Therefore, the same value as for
> :yref:`Bo1_Sphere_Aabb.aabbEnlargeFactor` must be given to it
> (:yref:`Ig2_Sphere_Sphere_ScGeom.interactionDetectionFactor` ).
>
> Note that only :yref:`Sphere` + :yref:`Sphere` interactions are
> supported; there is no parameter analogous to
> :yref:`distFactor<Ig2_Sphere_Sphere_ScGeom.interactionDetectionFactor>` in
> :yref:`Ig2_Facet_Sphere_ScGeom`. This is on purpose, since the interaction
> radius is meaningful in bulk material represented by sphere packing,
> whereas facets usually represent boundary conditions which should be exempt
> from this dense interaction network.
>
> @@ -576,12 +576,7 @@
>
> Again, missing combination will cause given shape combinations to
> freely interpenetrate one another.
>
> -#. :yref:`IGeom` type accepted by the ``Law2`` functor (below); it is the
> first part of functor's name after ``Law2`` (for instance,
> :yref:`Law2_ScGeom_CpmPhys_Cpm` accepts :yref:`ScGeom`). This is (for most
> cases) either :yref:`Dem3DofGeom` (total shear formulation) or
> :yref:`ScGeom` (incremental shear formulation). For :yref:`Dem3DofGeom`,
> the above example would simply change to::
> -
> - InteractionLoop([
> - Ig2_Sphere_Sphere_Dem3DofGeom(),
> - Ig2_Facet_Sphere_Dem3DofGeom()
> - ],[],[])
> +#. :yref:`IGeom` type accepted by the ``Law2`` functor (below); it is the
> first part of functor's name after ``Law2`` (for instance,
> :yref:`Law2_ScGeom_CpmPhys_Cpm` accepts :yref:`ScGeom`).
>
> Ip2 functors
> ^^^^^^^^^^^^
> @@ -693,29 +688,35 @@
> Imposing motion and forces
> --------------------------
>
> -* If a degree of freedom is blocked and a velocity is assigned along that
> direction (translational or rotational velocity), then the body will move
> at constant velocity. This is the simpler and recommended method to impose
> the motion of a body. This, for instance, will result in a constant
> velocity along $x$::
> +Imposed velocity
> +^^^^^^^^^^^^^^^^
> +
> +If a degree of freedom is blocked and a velocity is assigned along that
> direction (translational or rotational velocity), then the body will move
> at constant velocity. This is the simpler and recommended method to impose
> the motion of a body. This, for instance, will result in a constant
> velocity along $x$::
>
> O.bodies.append(utils.sphere((0,0,0),1))
> O.bodies[0].state.blockedDOFs='x'
> O.bodies[0].state.vel=(10,0,0)
>
> - Conversely, modifying the position directly is likely to break Yade's
> algorithms, especially those related to collision detection and contact
> laws, as they are based oon bodies velocities. Therefore, unless you really
> know what you are doing, don't do that for imposing a motion::
> +Conversely, modifying the position directly is likely to break Yade's
> algorithms, especially those related to collision detection and contact
> laws, as they are based oon bodies velocities. Therefore, unless you really
> know what you are doing, don't do that for imposing a motion::
>
> O.bodies.append(utils.sphere((0,0,0),1))
> O.bodies[0].state.blockedDOFs='x'
> O.bodies[0].state.pos=10*O.dt #REALLY BAD! Don't assign position
>
> -* Applying a force or a torque on a body is done via functions of the
> :yref:`ForceContainer`. It is as simple as this::
> +Imposed force
> +^^^^^^^^^^^^^
> +
> +Applying a force or a torque on a body is done via functions of the
> :yref:`ForceContainer`. It is as simple as this::
>
> O.forces.addF(0,(1,0,0)) #applies for one step
>
> - By default, the force applies for one time step only, and is resetted
> at the beginning of each step. For this reason, imposing a force at the
> begining of one step will have no effect at all, since it will be
> immediatly resetted. The only way is to place a :yref:`PyRunner` inside the
> simulation loop.
> +By default, the force applies for one time step only, and is resetted at
> the beginning of each step. For this reason, imposing a force at the
> begining of one step will have no effect at all, since it will be
> immediatly resetted. The only way is to place a :yref:`PyRunner` inside the
> simulation loop.
>
> - Applying the force permanently is possible with an optional argument
> (in this case it does not matter if the command comes at the begining of
> the time step)::
> +Applying the force permanently is possible with an optional argument (in
> this case it does not matter if the command comes at the begining of the
> time step)::
>
> O.forces.addF(0,(1,0,0),permanent=True) #applies permanently
>
> - The force will persist across iterations, until it is overwritten by
> another call to ``O.forces.addF(id,f,True)`` or erased by
> ``O.forces.reset(resetAll=True)``. The permanent force on a body can be
> checked with ``O.forces.permF(id)``.
> +The force will persist across iterations, until it is overwritten by
> another call to ``O.forces.addF(id,f,True)`` or erased by
> ``O.forces.reset(resetAll=True)``. The permanent force on a body can be
> checked with ``O.forces.permF(id)``.
>
> Boundary controllers
> --------------------
> @@ -1049,7 +1050,7 @@
> By editing the generated .gnuplot file you can plot any of the added Data
> afterwards.
>
>
> -* Data file is saved (compressed using bzip2) separately from the gnuplot
> file, so any other programs can be used to process them. In particular, the
> ``numpy.genfromtxt`` (documented `here <
> http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_)
> can be useful to import those data back to python; the decompression
> happens automatically.
> +* Data file is saved (compressed using bzip2) separately from the gnuplot
> file, so any other programs can be used to process them. In particular, the
> ``numpy.genfromtxt`` (`documented here <
> http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html>`_)
> can be useful to import those data back to python; the decompression
> happens automatically.
>
> * The gnuplot file can be run through gnuplot to produce the figure; see
> :yref:`yade.plot.saveGnuplot` documentation for details.
>
> @@ -1470,7 +1471,7 @@
> .. figure:: fig/micro-domains.*
>
> Micro-strain
> -""""""""""""
> +------------
> Below is an output of the :yref:`defToVtk<TesselationWrapper::defToVtk>`
> function visualized with paraview (in this case Yade's TesselationWrapper
> was used to process experimental data obtained on sand by Edward Ando at
> Grenoble University, 3SR lab.). The output is visualized with paraview, as
> explained in the previous section. Similar results can be generated from
> simulations:
>
> .. code-block:: python
> @@ -1493,7 +1494,7 @@
> .. figure:: fig/localstrain.*
>
> Micro-stress
> -""""""""""""
> +------------
> Stress fields can be generated by combining the volume returned by
> TesselationWrapper to per-particle stress given by
> :yref:`bodyStressTensors`. Since the stress $\sigma$ from bodyStressTensor
> implies a division by the volume $V_b$ of the solid particle, one has to
> re-normalize it in order to obtain the micro-stress as defined in
> [Catalano2013a]_ (equation 39 therein), i.e. $\overline{\sigma}^k =
> \sigma^k \times V_b^k / V_{\sigma}^k$ where $V_{\sigma}^k$ is the volume
> assigned to particle $k$ in the tesselation. For instance:
>
> .. code-block:: python
>
> === modified file 'doc/yade-articles.bib'
> --- doc/yade-articles.bib 2013-08-14 08:01:55 +0000
> +++ doc/yade-articles.bib 2013-08-26 21:57:21 +0000
> @@ -30,7 +30,7 @@
> url = {http://arxiv.org/pdf/1304.4895.pdf}
> }
>
> -@article {Chareyre2012a,
> +@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},
> @@ -359,15 +359,15 @@
> }
>
> @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}
> + 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,
>
> === modified file 'doc/yade-conferences.bib'
> --- doc/yade-conferences.bib 2013-08-14 08:01:55 +0000
> +++ doc/yade-conferences.bib 2013-08-26 21:53:51 +0000
> @@ -563,7 +563,7 @@
> doi = {10.1063/1.4812158}
> }
>
> -@InProceedings{ Hadda2013,
> +@InProceedings{ Hadda2013b,
> 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},
> @@ -590,3 +590,14 @@
> url = {https://www.yade-dem.org/publi/APC000495.pdf},
> doi = {10.1063/1.4811976}
> }
> +
> +@InProceedings{Maurin2013,
> + title={Discrete element modelling of bedload transport},
> + author={Maurin, Raphael and Chareyre, Bruno and Chauchat, Julien and
> Frey, Philippe},
> + booktitle={Proceedings of THESIS 2013, Two-pHase modElling for Sediment
> dynamIcS in Geophysical Flows},
> + adress = {Chatou, France},
> + pages={6p},
> + url = {https://www.yade-dem.org/publi/Maurin_et_al_THESIS_meta.pdf},
> + year={2013}
> +}
> +
>
> === modified file 'doc/yade-theses.bib'
> --- doc/yade-theses.bib 2013-08-19 17:51:24 +0000
> +++ doc/yade-theses.bib 2013-08-28 10:38:02 +0000
> @@ -101,6 +101,14 @@
> author = {Benoit Charlas},
> school = {Université de Grenoble},
> title = {Etude du comportement mécanique d'un hydrure
> intermétallique utilisé pour le stockage d'hydrogène},
> - year = {2013}
> + year = {2013},
> + url = {
> https://www.yade-dem.org/w/images/8/89/These_BenoitCharlas.pdf}
> }
>
> +@phdthesis{Catalano2012,
> + author = {Emanuele Catalano},
> + school = {Université de Grenoble},
> + title = {A pore-scale coupled hydromechanical model for biphasic
> granular media},
> + year = {2012},
> + url = {https://yade-dem.org/publi/Catalano_Thesis.pdf}
> +}
> \ No newline at end of file
>
> === removed file 'examples/test/CundallStrackTest.py'
> --- examples/test/CundallStrackTest.py 2013-03-28 11:13:12 +0000
> +++ examples/test/CundallStrackTest.py 1970-01-01 00:00:00 +0000
> @@ -1,21 +0,0 @@
> -# -*- coding: utf-8 -*-
> -
> -
> -O.engines=[
> - ForceResetter(),
> - InsertionSortCollider([Bo1_Sphere_Aabb(),]),
> - IGeomDispatcher([Ig2_Sphere_Sphere_Dem3DofGeom()]),
> - IPhysDispatcher([Ip2_2xFrictMat_CSPhys()]),
> - LawDispatcher([Law2_Dem3Dof_CSPhys_CundallStrack()]),
> - NewtonIntegrator(damping = 0.01,gravity=[0,0,-9.81])
> -]
> -
> -
> -
> -O.bodies.append(sphere([0,0,6],1,fixed=False, color=[0,1,0]))
> -O.bodies.append(sphere([0,0,0],1,fixed=True, color=[0,0,1]))
> -O.dt=.2*PWaveTimeStep()
> -
> -from yade import qt
> -qt.Controller()
> -qt.View()
>
> === removed file 'examples/test/Dem3DofGeom.py'
> --- examples/test/Dem3DofGeom.py 2013-03-28 16:15:54 +0000
> +++ examples/test/Dem3DofGeom.py 1970-01-01 00:00:00 +0000
> @@ -1,40 +0,0 @@
> -"Script showing shear interaction between facet/wall and sphere."
> -O.bodies.append([
> - #sphere([0,0,0],1,dynamic=False,color=(0,1,0),wire=True),
> -
> facet(([2,2,1],[-2,0,1],[2,-2,1]),fixed=True,color=(0,1,0),wire=False),
> - #wall([0,0,1],axis=2,color=(0,1,0)),
> - sphere([-1,0,2],1,fixed=False,color=(1,0,0),wire=True),
> -])
> -O.engines=[
> - ForceResetter(),
> - InsertionSortCollider([
> - Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()
> - ]),
> - IGeomDispatcher([
> - Ig2_Sphere_Sphere_Dem3DofGeom(),
> - Ig2_Facet_Sphere_Dem3DofGeom(),
> - Ig2_Wall_Sphere_Dem3DofGeom()
> - ]),
> - RotationEngine(rotationAxis=[0,1,0],angularVelocity=10,ids=[1]),
> - TranslationEngine(translationAxis=[1,0,0],velocity=10,ids=[1]),
> - NewtonIntegrator()#gravity=(0,0,-10))
> -]
> -O.miscParams=[
> -
> Gl1_Dem3DofGeom_SphereSphere(normal=True,rolledPoints=True,unrolledPoints=True,shear=True,shearLabel=True),
> -
> Gl1_Dem3DofGeom_FacetSphere(normal=False,rolledPoints=True,unrolledPoints=True,shear=True,shearLabel=True),
> -
> Gl1_Dem3DofGeom_WallSphere(normal=False,rolledPoints=True,unrolledPoints=True,shear=True,shearLabel=True),
> - Gl1_Sphere(wire=True)
> -]
> -
> -try:
> - from yade import qt
> - renderer=qt.Renderer()
> - renderer.wire=True
> - renderer.intrGeom=True
> - qt.Controller()
> - qt.View()
> -except ImportError: pass
> -
> -O.dt=1e-6
> -O.saveTmp('init')
> -O.run(1)
>
> === modified file 'pkg/common/Cylinder.hpp'
> --- pkg/common/Cylinder.hpp 2013-05-29 09:48:51 +0000
> +++ pkg/common/Cylinder.hpp 2013-08-29 10:30:31 +0000
> @@ -265,7 +265,7 @@
> public:
> //OpenMPAccumulator<Real> plasticDissipation;
> virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys,
> Interaction* I);
> -
> YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law
> for linear compression, and Mohr-Coulomb plasticity surface without
> cohesion.\nThis law implements the classical linear elastic-plastic law
> from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is
> (with the convention of positive tensile forces) $F_n=\\min(k_n u_n, 0)$.
> The shear force is $F_s=k_s u_s$, the plasticity condition defines the
> maximum value of the shear force : $F_s^{\\max}=F_n\\tan(\\phi)$, with
> $\\phi$ the friction angle.\n\n.. note::\n This law uses :yref:`ScGeom`;
> there is also functionally equivalent
> :yref:`Law2_Dem3DofGeom_FrictPhys_CundallStrack`, which uses
> :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the
> latest).\n\n.. note::\n This law is well tested in the context of triaxial
> simulation, and has been used for a number of published results (see e.g.
> [Scholtes2009b]_ and other papers from the same authors). It is generalised
> by :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`, which adds cohesion
> and moments at contact.",
> +
> YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law
> for linear compression, and Mohr-Coulomb plasticity surface without
> cohesion.\nThis law implements the classical linear elastic-plastic law
> from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is
> (with the convention of positive tensile forces) $F_n=\\min(k_n u_n, 0)$.
> The shear force is $F_s=k_s u_s$, the plasticity condition defines the
> maximum value of the shear force : $F_s^{\\max}=F_n\\tan(\\phi)$, with
> $\\phi$ the friction angle.\n\n.. note::\n This law is well tested in the
> context of triaxial simulation, and has been used for a number of published
> results (see e.g. [Scholtes2009b]_ and other papers from the same authors).
> It is generalised by :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`,
> which adds cohesion and moments at contact.",
> ((bool,neverErase,false,,"Keep
> interactions even if particles go away from each other (only in case
> another constitutive law is in the scene, e.g.
> :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
>
> ((bool,traceEnergy,false,Attr::hidden,"Define the total energy dissipated
> in plastic slips at all contacts."))
>
> ((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic
> dissipation (with O.trackEnergy)"))
>
> === modified file 'pkg/dem/ConcretePM.cpp'
> --- pkg/dem/ConcretePM.cpp 2013-08-23 15:21:20 +0000
> +++ pkg/dem/ConcretePM.cpp 2013-08-26 21:42:20 +0000
> @@ -345,7 +345,7 @@
> Vector3r& Fs(phys->Fs); /* for python access */
> const bool& isCohesive(phys->isCohesive);
>
> - Vector3r& epsTPl(phys->epsTPl);
> +// Vector3r& epsTPl(phys->epsTPl);
>
> #ifdef CPM_MATERIAL_MODEL
> Real& epsNPl(phys->epsNPl);
> @@ -589,7 +589,7 @@
> }
> }
> }
> - Real tr;
> +// Real tr;
> FOREACH(shared_ptr<Body> B, *scene->bodies){
> if (!B) continue;
> const Body::id_t& id = B->getId();
>
> === removed file 'pkg/dem/Dem3DofGeom_FacetSphere.cpp'
> --- pkg/dem/Dem3DofGeom_FacetSphere.cpp 2011-02-14 08:05:09 +0000
> +++ pkg/dem/Dem3DofGeom_FacetSphere.cpp 1970-01-01 00:00:00 +0000
> @@ -1,221 +0,0 @@
> -// © Václav Šmilauer <eudoxos@xxxxxxxx>
> -#include "Dem3DofGeom_FacetSphere.hpp"
> -#include<yade/pkg/common/Sphere.hpp>
> -#include<yade/pkg/common/Facet.hpp>
> -YADE_PLUGIN((Dem3DofGeom_FacetSphere)
> - #ifdef YADE_OPENGL
> - (Gl1_Dem3DofGeom_FacetSphere)
> - #endif
> - (Ig2_Facet_Sphere_Dem3DofGeom));
> -
> -CREATE_LOGGER(Dem3DofGeom_FacetSphere);
> -Dem3DofGeom_FacetSphere::~Dem3DofGeom_FacetSphere(){}
> -
> -void Dem3DofGeom_FacetSphere::setTgPlanePts(const Vector3r& p1new, const
> Vector3r& p2new){
> - TRVAR3(cp1pt,cp2rel,contPtInTgPlane2()-contPtInTgPlane1());
> -
> cp1pt=se31.orientation.conjugate()*(turnPlanePt(p1new,normal,se31.orientation*localFacetNormal)+contactPoint-se31.position);
> -
> cp2rel=se32.orientation.conjugate()*Dem3DofGeom_SphereSphere::rollPlanePtToSphere(p2new,effR2,-normal);
> - TRVAR3(cp1pt,cp2rel,contPtInTgPlane2()-contPtInTgPlane1());
> -}
> -
> -void Dem3DofGeom_FacetSphere::relocateContactPoints(const Vector3r& p1,
> const Vector3r& p2){
> - //TRVAR2(p2.norm(),effR2);
> - if(p2.squaredNorm()>pow(effR2,2)){
> - setTgPlanePts(Vector3r::Zero(),p2-p1);
> - }
> -}
> -
> -Real Dem3DofGeom_FacetSphere::slipToDisplacementTMax(Real
> displacementTMax){
> - //FIXME: not yet tested
> - // negative or zero: reset shear
> - if(displacementTMax<=0.){
> setTgPlanePts(Vector3r(0,0,0),Vector3r(0,0,0)); return displacementTMax;}
> - // otherwise
> - Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
> - Real currDistSq=(p2-p1).squaredNorm();
> - if(currDistSq<pow(displacementTMax,2)) return 0; // close enough,
> no slip needed
> - //Vector3r diff=.5*(sqrt(currDistSq)/displacementTMax-1)*(p2-p1);
> setTgPlanePts(p1+diff,p2-diff);
> - Real scale=displacementTMax/sqrt(currDistSq);
> setTgPlanePts(scale*p1,scale*p2);
> - return (displacementTMax/scale)*(1-scale);
> -}
> -
> -CREATE_LOGGER(Ig2_Facet_Sphere_Dem3DofGeom);
> -bool Ig2_Facet_Sphere_Dem3DofGeom::go(const shared_ptr<Shape>& cm1, const
> shared_ptr<Shape>& cm2, const State& state1, const State& state2, const
> Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
> - Facet* facet=static_cast<Facet*>(cm1.get());
> - Real sphereRadius=static_cast<Sphere*>(cm2.get())->radius;
> -
> - // IGeomFunctor::go(cm1,cm2,state1,state2,shift2,force,c);
> -
> - #if 1
> - /* new code written from scratch, to make sure the
> algorithm is correct; it is about the same speed
> - as sega's algo below, but seems more readable to
> me.
> - The FACET_TOPO thing is still missing here but can
> be copied literally once it is tested */
> - // begin facet-local coordinates
> - Vector3r
> cogLine=state1.ori.conjugate()*(state2.pos+shift2-state1.pos); // connect
> centers of gravity
> - //TRVAR4(state1.pos,state1.ori,state2.pos,cogLine);
> - Vector3r normal=facet->normal;
> - Real planeDist=normal.dot(cogLine);
> - if(planeDist<0){normal*=-1; planeDist*=-1; }
> - if(planeDist>sphereRadius && !c->isReal() &&
> !force) { /* LOG_TRACE("Sphere too far ("<<planeDist<<") from plane"); */
> return false; }
> - Vector3r planarPt=cogLine-planeDist*normal; //
> project sphere center to the facet plane
> - Real normDotPt[3];
> - Vector3r contactPt(Vector3r::Zero());
> - for(int i=0; i<3; i++)
> normDotPt[i]=facet->ne[i].dot(planarPt-facet->vertices[i]);
> - short
> w=(normDotPt[0]>0?1:0)+(normDotPt[1]>0?2:0)+(normDotPt[2]>0?4:0);
> -
> //TRVAR4(planarPt,normDotPt[0],normDotPt[1],normDotPt[2]);
> - //TRVAR2(normal,cogLine);
> -
> //TRVAR3(facet->vertices[0],facet->vertices[1],facet->vertices[2]);
> - switch(w){
> - case 0: contactPt=planarPt; break; //
> inside triangle
> - case 1:
> contactPt=getClosestSegmentPt(planarPt,facet->vertices[0],facet->vertices[1]);
> break; // +-- (n1)
> - case 2:
> contactPt=getClosestSegmentPt(planarPt,facet->vertices[1],facet->vertices[2]);
> break; // -+- (n2)
> - case 4:
> contactPt=getClosestSegmentPt(planarPt,facet->vertices[2],facet->vertices[0]);
> break; // --+ (n3)
> - case 3: contactPt=facet->vertices[1];
> break; // ++- (v1)
> - case 5: contactPt=facet->vertices[0];
> break; // +-+ (v0)
> - case 6: contactPt=facet->vertices[2];
> break; // -++ (v2)
> - case 7: throw logic_error("Impossible
> triangle intersection?"); // +++ (impossible)
> - default: throw logic_error("Nonsense
> intersection value!");
> - }
> - normal=cogLine-contactPt; // called normal, but it
> is no longer the facet's normal (for compat)
> - //TRVAR3(normal,contactPt,sphereRadius);
> - if(!c->isReal() &&
> normal.squaredNorm()>sphereRadius*sphereRadius && !force) { /*
> LOG_TRACE("Sphere too far from closest point"); */ return false; } // fast
> test before sqrt
> - Real norm=normal.norm(); normal/=norm; // normal
> is unit vector now
> - Real penetrationDepth=sphereRadius-norm;
> - #else
> - /* This code was mostly copied from
> InteractingFacet2InteractinSphere4SpheresContactGeometry */
> - // begin facet-local coordinates
> - Vector3r
> contactLine=state1.ori.Conjugate()*(state2.pos+shift2-state1.pos);
> - Vector3r normal=facet->normal;
> - Real L=normal.Dot(contactLine); // height/depth of
> sphere's center from facet's plane
> - if(L<0){normal*=-1; L*=-1;}
> - if(L>sphereRadius && !c->isReal()) return false;
> // sphere too far away from the plane
> -
> - Vector3r contactPt=contactLine-L*normal; //
> projection of sphere's center to facet's plane (preliminary contact point)
> - const Vector3r* edgeNormals=facet->ne; // array[3]
> of edge normals (in facet plane)
> - int edgeMax=0; Real
> distMax=edgeNormals[0].Dot(contactPt);
> - for(int i=1; i<3; i++){
> - Real dist=edgeNormals[i].Dot(contactPt);
> - if(distMax<dist){edgeMax=i; distMax=dist;}
> - }
> - //TRVAR2(distMax,edgeMax);
> - // OK, what's the logic here? Copying from
> IF2IS4SCG…
> - Real sphereRReduced=shrinkFactor*sphereRadius;
> - Real inCircleR=facet->icr-sphereRReduced;
> - Real penetrationDepth;
> - if(inCircleR<0){inCircleR=facet->icr;
> sphereRReduced=0;}
> - if(distMax<inCircleR){// contact with facet's
> surface
> - penetrationDepth=sphereRadius-L;
> - normal.normalize();
> - } else { // contact with the edge
> -
> contactPt+=edgeNormals[edgeMax]*(inCircleR-distMax);
> - bool noVertexContact=false;
> -
> //TRVAR3(edgeNormals[edgeMax],inCircleR,distMax);
> - // contact with vertex no. edgeMax
> - // FIXME: this is the original version,
> but why (edgeMax-1)%3? IN that case, edgeNormal to edgeMax would never be
> tried
> - // if (contactPt.Dot(edgeNormals[
> (edgeMax-1)%3])>inCircleR)
> contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
> - if (contactPt.Dot(edgeNormals[
> edgeMax ])>inCircleR)
> contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
> - // contact with vertex no. edgeMax+1
> - else
> if(contactPt.Dot(edgeNormals[edgeMax=(edgeMax+1)%3])>inCircleR)
> contactPt=facet->vu[edgeMax]*(facet->vl[edgeMax]-sphereRReduced);
> - // contact with edge no. edgeMax
> - else noVertexContact=true;
> - normal=contactLine-contactPt;
> - #ifdef FACET_TOPO
> - if(noVertexContact &&
> facet->edgeAdjIds[edgeMax]!=Body::ID_NONE){
> - // find angle between our
> normal and the facet's normal (still local coords)
> - Quaternionr q;
> q.Align(facet->normal,normal); AngleAxisr aa(q);
> - assert(aa.angle()>=0 &&
> aa.angle()<=Mathr::PI);
> -
> if(edgeNormals[edgeMax].Dot(aa.axis())<0) aa.angle()*=-1.;
> - bool
> negFace=normal.Dot(facet->normal)<0; // contact in on the negative facet's
> face
> - Real
> halfAngle=(negFace?-1.:1.)*facet->edgeAdjHalfAngle[edgeMax];
> - if(halfAngle<0 &&
> aa.angle()>halfAngle) return false; // on concave boundary, and if in the
> other facet's sector, no contact
> - // otherwise the contact
> will be created
> - }
> - #endif
> -
> //TRVAR4(contactLine,contactPt,normal,normal.norm());
> -
> //TRVAR3(se31.orientation*contactLine,se31.position+se31.orientation*contactPt,se31.orientation*normal);
> - Real norm=normal.norm(); normal/=norm;
> - penetrationDepth=sphereRadius-norm;
> - //TRVAR1(penetrationDepth);
> - }
> - // end facet-local coordinates
> - #endif
> -
> - if(penetrationDepth<0 && !c->isReal()) return false;
> -
> -
> - shared_ptr<Dem3DofGeom_FacetSphere> fs;
> - Vector3r normalGlob=state1.ori*normal;
> - bool isNew=false;
> - if(c->geom) fs=YADE_PTR_CAST<Dem3DofGeom_FacetSphere>(c->geom);
> - else {
> - // LOG_TRACE("Creating new Dem3DofGeom_FacetSphere");
> - fs=shared_ptr<Dem3DofGeom_FacetSphere>(new
> Dem3DofGeom_FacetSphere());
> - c->geom=fs;
> - isNew=true;
> - fs->effR2=sphereRadius-penetrationDepth;
> - fs->refR1=-1; fs->refR2=sphereRadius;
> - // postponed till below, to avoid numeric issues
> - // see https://lists.launchpad.net/yade-dev/msg02794.html
> - // since displacementN() is computed from fs->contactPoint,
> - // it was returning something like +1e-16 at the very
> first step
> - // when it was created ⇒ the constitutive law was erasing
> the
> - // contact as soon as it was created.
> - // fs->refLength=…
> - fs->cp1pt=contactPt; // facet-local intial contact point
> - fs->localFacetNormal=facet->normal;
> -
> fs->cp2rel.setFromTwoVectors(Vector3r::UnitX(),state2.ori.conjugate()*(-normalGlob));
> // initial sphere-local center-contactPt orientation WRT +x
> - fs->cp2rel.normalize();
> - }
> - fs->se31=state1.se3; fs->se32=state2.se3;
> fs->se32.position+=shift2;
> - fs->normal=normalGlob;
> -
> fs->contactPoint=state2.pos+shift2+(-normalGlob)*(sphereRadius-penetrationDepth);
> - // this refLength computation mimics what displacementN() does
> inside
> - // displcementN will therefore return exactly zero at the step the
> contact
> - // was created, which is what we want
> - if(isNew)
> fs->refLength=(state2.pos+shift2-fs->contactPoint).norm();
> - return true;
> -}
> -
> -#ifdef YADE_OPENGL
> -
> - #include<yade/lib/opengl/OpenGLWrapper.hpp>
> - #include<yade/lib/opengl/GLUtils.hpp>
> -
> - bool Gl1_Dem3DofGeom_FacetSphere::normal=false;
> - bool Gl1_Dem3DofGeom_FacetSphere::rolledPoints=false;
> - bool Gl1_Dem3DofGeom_FacetSphere::unrolledPoints=false;
> - bool Gl1_Dem3DofGeom_FacetSphere::shear=false;
> - bool Gl1_Dem3DofGeom_FacetSphere::shearLabel=false;
> -
> - void Gl1_Dem3DofGeom_FacetSphere::go(const shared_ptr<IGeom>& ig,
> const shared_ptr<Interaction>& ip, const shared_ptr<Body>& b1, const
> shared_ptr<Body>& b2, bool wireFrame){
> - Dem3DofGeom_FacetSphere* fs =
> static_cast<Dem3DofGeom_FacetSphere*>(ig.get());
> - const Se3r& se31=b1->state->se3,se32=b2->state->se3;
> - const Vector3r& pos1=se31.position; const Vector3r&
> pos2=se32.position;
> - const Quaternionr& ori1=se31.orientation; const
> Quaternionr& ori2=se32.orientation;
> - const Vector3r& contPt=fs->contactPoint;
> -
> - if(normal){
> -
> GLUtils::GLDrawArrow(contPt,contPt+fs->refLength*fs->normal); // normal of
> the contact
> - }
> - // sphere center to point on the sphere
> - if(rolledPoints){
> - //cerr<<pos1<<" "<<pos1+ori1*fs->cp1pt<<"
> "<<contPt<<endl;
> -
> GLUtils::GLDrawLine(pos1+ori1*fs->cp1pt,contPt,Vector3r(0,.5,1));
> -
> GLUtils::GLDrawLine(pos2,pos2+(ori2*fs->cp2rel*Vector3r::UnitX()*fs->effR2),Vector3r(0,1,.5));
> - }
> - // contact point to projected points
> - if(unrolledPoints||shear){
> - Vector3r ptTg1=fs->contPtInTgPlane1(),
> ptTg2=fs->contPtInTgPlane2();
> - if(unrolledPoints){
> - //TRVAR3(ptTg1,ptTg2,ss->normal)
> -
> GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1));
> -
> GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5));
> GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
> - }
> - if(shear){
> -
> GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
> - if(shearLabel)
> GLUtils::GLDrawNum(fs->displacementT().norm(),contPt,Vector3r(1,1,1));
> - }
> - }
> - }
> -
> -#endif
> -
>
> === removed file 'pkg/dem/Dem3DofGeom_FacetSphere.hpp'
> --- pkg/dem/Dem3DofGeom_FacetSphere.hpp 2013-08-23 15:21:20 +0000
> +++ pkg/dem/Dem3DofGeom_FacetSphere.hpp 1970-01-01 00:00:00 +0000
> @@ -1,80 +0,0 @@
> -#pragma once
> -#include<yade/pkg/dem/DemXDofGeom.hpp>
> -// for static roll/unroll functions in Dem3DofGeom_SphereSphere
> -#include<yade/pkg/dem/Dem3DofGeom_SphereSphere.hpp>
> -#include<yade/pkg/common/Sphere.hpp>
> -#include<yade/pkg/common/Facet.hpp>
> -
> -class Dem3DofGeom_FacetSphere: public Dem3DofGeom{
> - //! turn planePt in plane with normal0 around line passing through
> origin to plane with normal1
> - static Vector3r turnPlanePt(const Vector3r& planePt, const
> Vector3r& normal0, const Vector3r& normal1){ Quaternionr q;
> q.setFromTwoVectors(normal0,normal1); return q*planePt; }
> -
> - Vector3r contPtInTgPlane1() const { return
> turnPlanePt(se31.position+se31.orientation*cp1pt-contactPoint,se31.orientation*localFacetNormal,normal);
> }
> - Vector3r contPtInTgPlane2() const { return
> Dem3DofGeom_SphereSphere::unrollSpherePtToPlane(se32.orientation*cp2rel,effR2,-normal);}
> -
> - public:
> - virtual ~Dem3DofGeom_FacetSphere();
> - /******* API ********/
> - virtual Real displacementN(){ return
> (se32.position-contactPoint).norm()-refLength;}
> - virtual Vector3r displacementT(){ relocateContactPoints();
> return contPtInTgPlane2()-contPtInTgPlane1(); }
> - virtual Real slipToDisplacementTMax(Real displacementTMax);
> - /***** end API ******/
> -
> - void setTgPlanePts(const Vector3r&, const Vector3r&);
> - void relocateContactPoints(){
> relocateContactPoints(contPtInTgPlane1(),contPtInTgPlane2()); }
> - void relocateContactPoints(const Vector3r& p1, const
> Vector3r& p2);
> -
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(Dem3DofGeom_FacetSphere,Dem3DofGeom,"Class
> representing facet+sphere in contact which computes 3 degrees of freedom
> (normal and shear deformation).",
> - ((Vector3r,cp1pt,,,"Reference contact point on the facet
> in facet-local coords."))
> - ((Quaternionr,cp2rel,,,"Orientation between +x and the
> reference contact point (on the sphere) in sphere-local coords"))
> - ((Vector3r,localFacetNormal,,,"Unit normal of the facet
> plane in facet-local coordinates"))
> - ((Real,effR2,,,"Effective radius of sphere")),
> - /*ctor*/ createIndex();
> - );
> - REGISTER_CLASS_INDEX(Dem3DofGeom_FacetSphere,Dem3DofGeom);
> - DECLARE_LOGGER;
> - friend class Gl1_Dem3DofGeom_FacetSphere;
> - friend class Ig2_Facet_Sphere_Dem3DofGeom;
> -};
> -REGISTER_SERIALIZABLE(Dem3DofGeom_FacetSphere);
> -
> -#ifdef YADE_OPENGL
> - #include<yade/pkg/common/GLDrawFunctors.hpp>
> - class Gl1_Dem3DofGeom_FacetSphere:public GlIGeomFunctor{
> - public:
> - virtual void go(const shared_ptr<IGeom>&,const
> shared_ptr<Interaction>&,const shared_ptr<Body>&,const
> shared_ptr<Body>&,bool wireFrame);
> - RENDERS(Dem3DofGeom_FacetSphere);
> -
> YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_Dem3DofGeom_FacetSphere,GlIGeomFunctor,"Render
> interaction of facet and sphere (represented by Dem3DofGeom_FacetSphere)",
> - ((bool,normal,false,,"Render interaction normal"))
> - ((bool,rolledPoints,false,,"Render points rolled
> on the sphere & facet (original contact point)"))
> - ((bool,unrolledPoints,false,,"Render original
> contact points unrolled to the contact plane"))
> - ((bool,shear,false,,"Render shear line in the
> contact plane"))
> - ((bool,shearLabel,false,,"Render shear magnitude
> as number"))
> - );
> - };
> - REGISTER_SERIALIZABLE(Gl1_Dem3DofGeom_FacetSphere);
> -#endif
> -
> -#include<yade/pkg/common/Dispatching.hpp>
> -class Ig2_Facet_Sphere_Dem3DofGeom:public IGeomFunctor{
> - Vector3r getClosestSegmentPt(const Vector3r& P, const Vector3r& A,
> const Vector3r& B){
> - // algo:
> http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/
> - Vector3r BA=B-A;
> - Real u=(P.dot(BA)-A.dot(BA))/(BA.squaredNorm());
> - return A+min((Real)1.,max((Real)0.,u))*BA;
> - }
> - public:
> - virtual bool go(const shared_ptr<Shape>& cm1, const
> shared_ptr<Shape>& cm2, const State& state1, const State& state2, const
> Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
> - virtual bool goReverse( const shared_ptr<Shape>& cm1,
> const shared_ptr<Shape>& cm2, const State& state1, const State& state2,
> const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>&
> c){
> - c->swapOrder(); return
> go(cm2,cm1,state2,state1,-shift2,force,c);
> - LOG_ERROR("!! goReverse maybe doesn't work in
> Ig2_Facet_Sphere_Dem3DofGeom. IGeomDispatcher should swap interaction
> members first and call go(...) afterwards.");
> - }
> -
> - FUNCTOR2D(Facet,Sphere);
> - DEFINE_FUNCTOR_ORDER_2D(Facet,Sphere);
> -
> YADE_CLASS_BASE_DOC_ATTRS(Ig2_Facet_Sphere_Dem3DofGeom,IGeomFunctor,"Compute
> geometry of facet-sphere contact with normal and shear DOFs. As in all
> other Dem3DofGeom-related classes, total formulation of both shear and
> normal deformations is used. See :yref:`Dem3DofGeom_FacetSphere` for more
> information.",
> - // unused: ((Real,shrinkFactor,0.,,"Reduce the facet's
> size, probably to avoid singularities at common facets' edges (?)"))
> - );
> - DECLARE_LOGGER;
> -};
> -REGISTER_SERIALIZABLE(Ig2_Facet_Sphere_Dem3DofGeom);
> -
>
> === removed file 'pkg/dem/Dem3DofGeom_SphereSphere.cpp'
> --- pkg/dem/Dem3DofGeom_SphereSphere.cpp 2010-12-13 12:11:43 +0000
> +++ pkg/dem/Dem3DofGeom_SphereSphere.cpp 1970-01-01 00:00:00 +0000
> @@ -1,196 +0,0 @@
> -#include "Dem3DofGeom_SphereSphere.hpp"
> -
> -#include<yade/pkg/common/Sphere.hpp>
> -#include<yade/core/Omega.hpp>
> -YADE_PLUGIN((Dem3DofGeom_SphereSphere)
> - #ifdef YADE_OPENGL
> - (Gl1_Dem3DofGeom_SphereSphere)
> - #endif
> - (Ig2_Sphere_Sphere_Dem3DofGeom));
> -
> -
> -Dem3DofGeom_SphereSphere::~Dem3DofGeom_SphereSphere(){}
> -
> -/*! Project point from sphere surface to tangent plane,
> - * such that the angle of shortest arc from (1,0,0) pt on the sphere to
> the point itself is the same
> - * as the angle of segment of the same length on the tangent plane.
> - *
> - * This function is (or should be) inverse of ScGeom::rollPlanePtToSphere.
> - *
> - * @param fromXtoPtOri gives orientation of the vector from sphere center
> to the sphere point from the global +x axis.
> - * @param radius the distance from sphere center to the contact plane
> - * @param planeNormal unit vector pointing away from the sphere center,
> determining plane orientation on which the projected point lies.
> - * @returns The projected point coordinates (with origin at the contact
> point).
> - */
> -Vector3r Dem3DofGeom_SphereSphere::unrollSpherePtToPlane(const
> Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& planeNormal){
> - Quaternionr normal2pt;
> normal2pt.setFromTwoVectors(planeNormal,fromXtoPtOri*Vector3r::UnitX());
> - AngleAxisr aa(normal2pt);
> - return (aa.angle()*radius) /* length */
> *(aa.axis().cross(planeNormal)) /* direction: both are unit vectors */;
> -}
> -
> -/*! Project point from tangent plane to the sphere.
> - *
> - * This function is (or should be) inverse of
> ScGeom::unrollSpherePtToPlane.
> - *
> - * @param planePt point on the tangent plane, with origin at the contact
> point (i.e. at sphere center + normal*radius)
> - * @param radius sphere radius
> - * @param planeNormal _unit_ vector pointing away from sphere center
> - * @returns orientation that transforms +x axis to the vector between
> sphere center and point on the sphere that corresponds to planePt.
> - *
> - * @note It is not checked whether planePt relly lies on the tangent
> plane. If not, result will be incorrect.
> - */
> -Quaternionr Dem3DofGeom_SphereSphere::rollPlanePtToSphere(const Vector3r&
> planePt, const Real& radius, const Vector3r& planeNormal){
> - if (planePt!=Vector3r::Zero()) {
> - Quaternionr normal2pt;
> - Vector3r axis=planeNormal.cross(planePt); axis.normalize();
> - Real angle=planePt.norm()/radius;
> - normal2pt=Quaternionr(AngleAxisr(angle,axis));
> - Quaternionr ret; return
> ret.setFromTwoVectors(Vector3r::UnitX(),normal2pt*planeNormal);
> - } else {
> - Quaternionr ret; return
> ret.setFromTwoVectors(Vector3r::UnitX(),planeNormal);
> - }
> -}
> -
> -
> -
> -/* Set contact points on both spheres such that their projection is the
> one given
> - * (should be on the plane passing through origin and oriented with
> normal; not checked!)
> - */
> -void Dem3DofGeom_SphereSphere::setTgPlanePts(Vector3r p1new, Vector3r
> p2new){
> - cp1rel=ori1.conjugate()*rollPlanePtToSphere(p1new,effR1,normal);
> - cp2rel=ori2.conjugate()*rollPlanePtToSphere(p2new,effR2,-normal);
> -}
> -
> -
> -
> -/*! Perform slip of the projected contact points so that their distance
> becomes equal (or remains smaller) than the given one.
> - * The slipped distance is returned.
> - */
> -Real Dem3DofGeom_SphereSphere::slipToDisplacementTMax(Real
> displacementTMax){
> - // very close, reset shear
> - if(displacementTMax<=0.){
> setTgPlanePts(Vector3r(0,0,0),Vector3r(0,0,0)); return displacementTMax;}
> - // otherwise
> - Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
> - Real currDistSq=(p2-p1).squaredNorm();
> - if(currDistSq<pow(displacementTMax,2)) return 0; // close enough,
> no slip needed
> - Vector3r diff=.5*(displacementTMax/sqrt(currDistSq)-1)*(p2-p1);
> - setTgPlanePts(p1-diff,p2+diff);
> - return 2*diff.norm();
> -}
> -
> -
> -/* Move contact point on both spheres in such way that their relative
> position (displacementT) is the same;
> - * this should be done regularly to ensure that the angle doesn't go over
> π, since then quaternion would
> - * flip axis and the point would project on other side of the tangent
> plane piece. */
> -void Dem3DofGeom_SphereSphere::relocateContactPoints(){
> - relocateContactPoints(contPtInTgPlane1(),contPtInTgPlane2());
> -}
> -
> -/*! Like Dem3DofGeom_SphereSphere::relocateContactPoints(), but use
> already computed tangent plane points. */
> -void Dem3DofGeom_SphereSphere::relocateContactPoints(const Vector3r& p1,
> const Vector3r& p2){
> - Vector3r midPt=(effR1/(effR1+effR2))*(p1+p2); // proportionally to
> radii, so that angle would be the same
> - if((p1.squaredNorm()>pow(effR1,2) ||
> p2.squaredNorm()>pow(effR2,2)) &&
> midPt.squaredNorm()>pow(min(effR1,effR2),2)){
> - //cerr<<"RELOCATION with displacementT="<<displacementT();
> // should be the same before and after relocation
> - setTgPlanePts(p1-midPt,p2-midPt);
> - //cerr<<" → "<<displacementT()<<endl;
> - }
> -}
> -
> -
> -
> -#ifdef YADE_OPENGL
> - #include<yade/lib/opengl/OpenGLWrapper.hpp>
> - #include<yade/lib/opengl/GLUtils.hpp>
> - bool Gl1_Dem3DofGeom_SphereSphere::normal=false;
> - bool Gl1_Dem3DofGeom_SphereSphere::rolledPoints=false;
> - bool Gl1_Dem3DofGeom_SphereSphere::unrolledPoints=false;
> - bool Gl1_Dem3DofGeom_SphereSphere::shear=false;
> - bool Gl1_Dem3DofGeom_SphereSphere::shearLabel=false;
> -
> - void Gl1_Dem3DofGeom_SphereSphere::go(const shared_ptr<IGeom>& ig,
> const shared_ptr<Interaction>& ip, const shared_ptr<Body>& b1, const
> shared_ptr<Body>& b2, bool wireFrame){
> - Dem3DofGeom_SphereSphere* ss =
> static_cast<Dem3DofGeom_SphereSphere*>(ig.get());
> - //const Se3r&
> se31=b1->physicalParameters->dispSe3,se32=b2->physicalParameters->dispSe3;
> - const Se3r& se31=b1->state->se3,se32=b2->state->se3;
> - const Vector3r& pos1=se31.position,pos2=se32.position;
> - Vector3r& contPt=ss->contactPoint;
> -
> - if(normal){
> -
> GLUtils::GLDrawArrow(contPt,contPt+ss->normal*.5*ss->refLength); // normal
> of the contact
> - }
> - #if 0
> - // never used, since bending/torsion not used
> - //Vector3r
> contPt=se31.position+(ss->effR1/ss->refLength)*(se32.position-se31.position);
> // must be recalculated to not be unscaled if scaling displacements ...
> - GLUtils::GLDrawLine(pos1,pos2,Vector3r(.5,.5,.5));
> - Vector3r bend; Real tors;
> - ss->bendingTorsionRel(bend,tors);
> -
> GLUtils::GLDrawLine(contPt,contPt+10*ss->radius1*(bend+ss->normal*tors),Vector3r(1,0,0));
> - #if 0
> -
> GLUtils::GLDrawNum(bend[0],contPt-.2*ss->normal*ss->radius1,Vector3r(1,0,0));
> -
> GLUtils::GLDrawNum(bend[1],contPt,Vector3r(0,1,0));
> -
> GLUtils::GLDrawNum(bend[2],contPt+.2*ss->normal*ss->radius1,Vector3r(0,0,1));
> -
> GLUtils::GLDrawNum(tors,contPt+.5*ss->normal*ss->radius2,Vector3r(1,1,0));
> - #endif
> - #endif
> - // sphere center to point on the sphere
> - if(rolledPoints){
> -
> GLUtils::GLDrawLine(pos1,pos1+(ss->ori1*ss->cp1rel*Vector3r::UnitX()*ss->effR1),Vector3r(0,.5,1));
> -
> GLUtils::GLDrawLine(pos2,pos2+(ss->ori2*ss->cp2rel*Vector3r::UnitX()*ss->effR2),Vector3r(0,1,.5));
> - }
> - //TRVAR4(pos1,ss->ori1,pos2,ss->ori2);
> -
> //TRVAR2(ss->cp2rel,pos2+(ss->ori2*ss->cp2rel*Vector3r::UnitX()*ss->effR2));
> - // contact point to projected points
> - if(unrolledPoints||shear){
> - Vector3r ptTg1=ss->contPtInTgPlane1(),
> ptTg2=ss->contPtInTgPlane2();
> - if(unrolledPoints){
> - //TRVAR3(ptTg1,ptTg2,ss->normal)
> -
> GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1));
> GLUtils::GLDrawLine(pos1,contPt+ptTg1,Vector3r(0,.5,1));
> -
> GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5));
> GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
> - }
> - if(shear){
> -
> GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
> - if(shearLabel)
> GLUtils::GLDrawNum(ss->displacementT().norm(),contPt,Vector3r(1,1,1));
> - }
> - }
> - }
> -#endif
> -
> -CREATE_LOGGER(Ig2_Sphere_Sphere_Dem3DofGeom);
> -
> -bool Ig2_Sphere_Sphere_Dem3DofGeom::go(const shared_ptr<Shape>& cm1,
> const shared_ptr<Shape>& cm2, const State& state1, const State& state2,
> const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>&
> c){
> - Sphere *s1=static_cast<Sphere*>(cm1.get()),
> *s2=static_cast<Sphere*>(cm2.get());
> - Vector3r normal=(state2.pos+shift2)-state1.pos;
> - Real
> penetrationDepthSq=pow((distFactor>0?distFactor:1.)*(s1->radius+s2->radius),2)-normal.squaredNorm();
> - if (penetrationDepthSq<0 && !c->isReal() && !force){
> - return false;
> - }
> -
> - Real dist=normal.norm(); normal/=dist; /* normal is unit vector
> now */
> - shared_ptr<Dem3DofGeom_SphereSphere> ss;
> - if(c->geom) ss=YADE_PTR_CAST<Dem3DofGeom_SphereSphere>(c->geom);
> - else {
> - ss=shared_ptr<Dem3DofGeom_SphereSphere>(new
> Dem3DofGeom_SphereSphere());
> - c->geom=ss;
> -
> - // constants
> - if(distFactor>0) ss->refLength=dist;
> - else ss->refLength=s1->radius+s2->radius;
> - ss->refR1=s1->radius; ss->refR2=s2->radius;
> -
> - Real penetrationDepth=s1->radius+s2->radius-ss->refLength;
> -
> - if(scene->iter<=10){
> - ss->effR1=s1->radius-.5*penetrationDepth;
> ss->effR2=s2->radius-.5*penetrationDepth;
> - } else {ss->effR1=s1->radius; ss->effR2=s2->radius;}
> -
> - // for bending only:
> ss->initRelOri12=state1.ori.Conjugate()*state2.ori;
> - // quasi-constants
> -
> ss->cp1rel.setFromTwoVectors(Vector3r::UnitX(),state1.ori.conjugate()*normal);
> -
> ss->cp2rel.setFromTwoVectors(Vector3r::UnitX(),state2.ori.conjugate()*(-normal));
> - ss->cp1rel.normalize(); ss->cp2rel.normalize();
> - }
> - ss->normal=normal;
> -
> ss->contactPoint=state1.pos+(ss->effR1-.5*(ss->refLength-dist))*ss->normal;
> - ss->se31=state1.se3; ss->se32=state2.se3;
> ss->se32.position+=shift2;
> - return true;
> -}
> -
>
> === removed file 'pkg/dem/Dem3DofGeom_SphereSphere.hpp'
> --- pkg/dem/Dem3DofGeom_SphereSphere.hpp 2011-01-09 16:34:50 +0000
> +++ pkg/dem/Dem3DofGeom_SphereSphere.hpp 1970-01-01 00:00:00 +0000
> @@ -1,80 +0,0 @@
> -#pragma once
> -#include<yade/pkg/dem/DemXDofGeom.hpp>
> -#include<yade/pkg/common/Sphere.hpp>
> -
> -class Dem3DofGeom_SphereSphere: public Dem3DofGeom{
> - public:
> - // auxiliary functions; the quaternion magic is all in here
> - static Vector3r unrollSpherePtToPlane(const Quaternionr&
> fromXtoPtOri, const Real& radius, const Vector3r& normal);
> - static Quaternionr rollPlanePtToSphere(const Vector3r&
> planePt, const Real& radius, const Vector3r& normal);
> - private:
> - Vector3r contPtInTgPlane1() const { return
> unrollSpherePtToPlane(ori1*cp1rel,effR1,normal); }
> - Vector3r contPtInTgPlane2() const { return
> unrollSpherePtToPlane(ori2*cp2rel,effR2,-normal);}
> - void setTgPlanePts(Vector3r p1new, Vector3r p2new);
> - void relocateContactPoints();
> - void relocateContactPoints(const Vector3r& tgPlanePt1,
> const Vector3r& tgPlanePt2);
> - public:
> - //! shorthands
> - const Vector3r &pos1; const Quaternionr &ori1; const
> Vector3r &pos2; const Quaternionr &ori2;
> - virtual ~Dem3DofGeom_SphereSphere();
> -
> - /********* API **********/
> - Real displacementN(){ return (pos2-pos1).norm()-refLength;
> }
> - Vector3r displacementT() {
> - // enabling automatic relocation decreases overall
> simulation speed by about 3%; perhaps: bool largeStrains ... ?
> - #if 1
> - Vector3r p1=contPtInTgPlane1(),
> p2=contPtInTgPlane2(); relocateContactPoints(p1,p2); return p2-p1; // shear
> before relocation, but that is OK
> - #else
> - return
> contPtInTgPlane2()-contPtInTgPlane1();
> - #endif
> - }
> - virtual Real slipToDisplacementTMax(Real displacementTMax);
> - /********* end API ***********/
> -
> -
> YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(Dem3DofGeom_SphereSphere,Dem3DofGeom,"Class
> representing 2 spheres in contact which computes 3 degrees of freedom
> (normal and shear deformation).",
> - ((Real,effR1,,,"Effective radius of sphere #1; can be
> smaller/larger than refR1 (the actual radius), but quasi-constant
> throughout interaction life"))
> - ((Real,effR2,,,"Same as effR1, but for sphere #2."))
> - ((Quaternionr,cp1rel,,,"Sphere's #1 relative orientation
> of the contact point with regards to sphere-local +x axis
> (quasi-constant)"))
> - ((Quaternionr,cp2rel,,,"Same as cp1rel, but for sphere
> #2.")),
> - /* extra initializers */
> ((pos1,se31.position))((ori1,se31.orientation))((pos2,se32.position))((ori2,se32.orientation)),
> - /*ctor*/ createIndex(); ,
> - /*py*/
> - );
> - REGISTER_CLASS_INDEX(Dem3DofGeom_SphereSphere,Dem3DofGeom);
> - friend class Gl1_Dem3DofGeom_SphereSphere;
> - friend class Ig2_Sphere_Sphere_Dem3DofGeom;
> -};
> -REGISTER_SERIALIZABLE(Dem3DofGeom_SphereSphere);
> -
> -#ifdef YADE_OPENGL
> - #include<yade/pkg/common/GLDrawFunctors.hpp>
> - class Gl1_Dem3DofGeom_SphereSphere:public GlIGeomFunctor{
> - public:
> - virtual void go(const shared_ptr<IGeom>&,const
> shared_ptr<Interaction>&,const shared_ptr<Body>&,const
> shared_ptr<Body>&,bool wireFrame);
> - RENDERS(Dem3DofGeom_SphereSphere);
> -
> YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_Dem3DofGeom_SphereSphere,GlIGeomFunctor,"Render
> interaction of 2 spheres (represented by Dem3DofGeom_SphereSphere)",
> - ((bool,normal,false,,"Render interaction normal"))
> - ((bool,rolledPoints,false,,"Render points rolled
> on the spheres (tracks the original contact point)"))
> - ((bool,unrolledPoints,false,,"Render original
> contact points unrolled to the contact plane"))
> - ((bool,shear,false,,"Render shear line in the
> contact plane"))
> - ((bool,shearLabel,false,,"Render shear magnitude
> as number"))
> - );
> - };
> - REGISTER_SERIALIZABLE(Gl1_Dem3DofGeom_SphereSphere);
> -#endif
> -
> -#include<yade/pkg/common/Dispatching.hpp>
> -class Ig2_Sphere_Sphere_Dem3DofGeom:public IGeomFunctor{
> - public:
> - virtual bool go(const shared_ptr<Shape>& cm1, const
> shared_ptr<Shape>& cm2, const State& state1, const State& state2, const
> Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
> - virtual bool goReverse( const shared_ptr<Shape>&, const
> shared_ptr<Shape>&, const State&, const State&, const Vector3r& shift2,
> const bool& force, const shared_ptr<Interaction>&){throw
> runtime_error("goReverse on symmetric functor should never be called!");}
> - FUNCTOR2D(Sphere,Sphere);
> - DEFINE_FUNCTOR_ORDER_2D(Sphere,Sphere);
> - DECLARE_LOGGER;
> -
> YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_Sphere_Dem3DofGeom,IGeomFunctor,
> - "Functor handling contact of 2 spheres, producing
> Dem3DofGeom instance",
> - ((Real,distFactor,-1,,"Factor of sphere radius such that
> sphere \"touch\" if their centers are not further than distFactor*(r1+r2);
> if negative, equilibrium distance is the sum of the sphere's radii."))
> - );
> -};
> -REGISTER_SERIALIZABLE(Ig2_Sphere_Sphere_Dem3DofGeom);
> -
>
> === removed file 'pkg/dem/Dem3DofGeom_WallSphere.cpp'
> --- pkg/dem/Dem3DofGeom_WallSphere.cpp 2010-12-13 12:11:43 +0000
> +++ pkg/dem/Dem3DofGeom_WallSphere.cpp 1970-01-01 00:00:00 +0000
> @@ -1,117 +0,0 @@
> -// © 2009 Václav Šmilauer <eudoxos@xxxxxxxx>
> -#include<yade/pkg/dem/Dem3DofGeom_WallSphere.hpp>
> -#include<yade/pkg/common/Sphere.hpp>
> -#include<yade/pkg/common/Wall.hpp>
> -YADE_PLUGIN((Dem3DofGeom_WallSphere)
> - #ifdef YADE_OPENGL
> - (Gl1_Dem3DofGeom_WallSphere)
> - #endif
> - (Ig2_Wall_Sphere_Dem3DofGeom));
> -
> -CREATE_LOGGER(Dem3DofGeom_WallSphere);
> -Dem3DofGeom_WallSphere::~Dem3DofGeom_WallSphere(){}
> -
> -void Dem3DofGeom_WallSphere::setTgPlanePts(const Vector3r& p1new, const
> Vector3r& p2new){
> - TRVAR3(cp1pt,cp2rel,contPtInTgPlane2()-contPtInTgPlane1());
> - cp1pt=p1new+contactPoint-se31.position;
> -
> cp2rel=se32.orientation.conjugate()*Dem3DofGeom_SphereSphere::rollPlanePtToSphere(p2new,effR2,-normal);
> - TRVAR3(cp1pt,cp2rel,contPtInTgPlane2()-contPtInTgPlane1());
> -}
> -
> -void Dem3DofGeom_WallSphere::relocateContactPoints(const Vector3r& p1,
> const Vector3r& p2){
> - //TRVAR2(p2.norm(),effR2);
> - if(p2.squaredNorm()>pow(effR2,2)){
> - setTgPlanePts(Vector3r::Zero(),p2-p1);
> - }
> -}
> -
> -Real Dem3DofGeom_WallSphere::slipToDisplacementTMax(Real
> displacementTMax){
> - //FIXME: not yet tested
> - // negative or zero: reset shear
> - if(displacementTMax<=0.){
> setTgPlanePts(Vector3r(0,0,0),Vector3r(0,0,0)); return displacementTMax;}
> - // otherwise
> - Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
> - Real currDistSq=(p2-p1).squaredNorm();
> - if(currDistSq<pow(displacementTMax,2)) return 0; // close enough,
> no slip needed
> - //Vector3r diff=.5*(sqrt(currDistSq)/displacementTMax-1)*(p2-p1);
> setTgPlanePts(p1+diff,p2-diff);
> - Real scale=displacementTMax/sqrt(currDistSq);
> setTgPlanePts(scale*p1,scale*p2);
> - return (displacementTMax/scale)*(1-scale);
> -}
> -
> -CREATE_LOGGER(Ig2_Wall_Sphere_Dem3DofGeom);
> -bool Ig2_Wall_Sphere_Dem3DofGeom::go(const shared_ptr<Shape>& cm1, const
> shared_ptr<Shape>& cm2, const State& state1, const State& state2, const
> Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
> - Wall* wall=static_cast<Wall*>(cm1.get());
> - Real sphereRadius=static_cast<Sphere*>(cm2.get())->radius;
> -
> - Real dist=(state2.pos+shift2)[wall->axis]-state1.pos[wall->axis];
> - if(!c->isReal() && abs(dist)>sphereRadius && !force){
> /*LOG_DEBUG("dist="<<dist<<", returning false");*/ return false; } // wall
> and sphere too far from each other
> -
> - // contact point is sphere center projected onto the wall
> - Vector3r contPt=state2.pos;
> contPt[wall->axis]=state1.pos[wall->axis];
> - Vector3r normalGlob(0.,0.,0.);
> - // wall interacting from both sides: normal depends on sphere's
> position
> - assert(wall->sense==-1 || wall->sense==0 || wall->sense==1);
> - if(wall->sense==0) normalGlob[wall->axis]=dist>0?1.:-1.;
> - else normalGlob[wall->axis]=wall->sense==1?1.:-1;
> -
> - shared_ptr<Dem3DofGeom_WallSphere> ws;
> - if(c->geom) ws=YADE_PTR_CAST<Dem3DofGeom_WallSphere>(c->geom);
> - else {
> - ws=shared_ptr<Dem3DofGeom_WallSphere>(new
> Dem3DofGeom_WallSphere());
> - c->geom=ws;
> - ws->effR2=abs(dist);
> - ws->refR1=-1; ws->refR2=sphereRadius;
> - ws->refLength=ws->effR2;
> - ws->cp1pt=contPt-state1.pos; // initial contact point
> relative to wall position (orientation is global, since it is coincident
> with local for a wall)
> - ws->cp2rel=Quaternionr::Identity();
> -
> ws->cp2rel.setFromTwoVectors(Vector3r::UnitX(),state2.ori.conjugate()*(-normalGlob));
> // initial sphere-local center-contactPt orientation WRT +x
> - ws->cp2rel.normalize();
> - //LOG_INFO(ws->cp1pt);
> - }
> - ws->se31=state1.se3; ws->se32=state2.se3;
> ws->se32.position+=shift2;
> - ws->contactPoint=contPt;
> - ws->normal=normalGlob;
> - return true;
> -}
> -#ifdef YADE_OPENGL
> -
> - #include<yade/lib/opengl/OpenGLWrapper.hpp>
> - #include<yade/lib/opengl/GLUtils.hpp>
> -
> - bool Gl1_Dem3DofGeom_WallSphere::normal=false;
> - bool Gl1_Dem3DofGeom_WallSphere::rolledPoints=false;
> - bool Gl1_Dem3DofGeom_WallSphere::unrolledPoints=false;
> - bool Gl1_Dem3DofGeom_WallSphere::shear=false;
> - bool Gl1_Dem3DofGeom_WallSphere::shearLabel=false;
> -
> - void Gl1_Dem3DofGeom_WallSphere::go(const shared_ptr<IGeom>& ig,
> const shared_ptr<Interaction>& ip, const shared_ptr<Body>& b1, const
> shared_ptr<Body>& b2, bool wireFrame){
> - Dem3DofGeom_WallSphere*
> ws=static_cast<Dem3DofGeom_WallSphere*>(ig.get());
> - const Se3r& se31=b1->state->se3,se32=b2->state->se3;
> - const Vector3r& pos1=se31.position; const Vector3r&
> pos2=se32.position;
> - //const Quaternionr& ori1=se31.orientation;
> - const Quaternionr& ori2=se32.orientation;
> - const Vector3r& contPt=ws->contactPoint;
> -
> - if(normal){
> -
> GLUtils::GLDrawArrow(contPt,contPt+ws->refLength*ws->normal); // normal of
> the contact
> - }
> - // sphere center to point on the sphere
> - if(rolledPoints){
> -
> GLUtils::GLDrawLine(pos1+ws->cp1pt,contPt,Vector3r(0,.5,1));
> -
> GLUtils::GLDrawLine(pos2,pos2+(ori2*ws->cp2rel*Vector3r::UnitX()*ws->effR2),Vector3r(0,1,.5));
> - }
> - // contact point to projected points
> - if(unrolledPoints||shear){
> - Vector3r ptTg1=ws->contPtInTgPlane1(),
> ptTg2=ws->contPtInTgPlane2();
> - if(unrolledPoints){
> -
> GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1));
> -
> GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5));
> GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
> - }
> - if(shear){
> -
> GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
> - if(shearLabel)
> GLUtils::GLDrawNum(ws->displacementT().norm(),contPt,Vector3r(1,1,1));
> - }
> - }
> - }
> -
> -#endif
>
> === removed file 'pkg/dem/Dem3DofGeom_WallSphere.hpp'
> --- pkg/dem/Dem3DofGeom_WallSphere.hpp 2011-01-09 16:34:50 +0000
> +++ pkg/dem/Dem3DofGeom_WallSphere.hpp 1970-01-01 00:00:00 +0000
> @@ -1,70 +0,0 @@
> -#pragma once
> -#include<yade/pkg/dem/DemXDofGeom.hpp>
> -// for static roll/unroll functions in Dem3DofGeom_SphereSphere
> -#include<yade/pkg/dem/Dem3DofGeom_SphereSphere.hpp>
> -#include<yade/pkg/common/Sphere.hpp>
> -#include<yade/pkg/common/Wall.hpp>
> -
> -class Dem3DofGeom_WallSphere: public Dem3DofGeom{
> -
> - Vector3r contPtInTgPlane1() const { return
> se31.position+cp1pt-contactPoint; }
> - Vector3r contPtInTgPlane2() const { return
> Dem3DofGeom_SphereSphere::unrollSpherePtToPlane(se32.orientation*cp2rel,effR2,-normal);}
> -
> - public:
> - virtual ~Dem3DofGeom_WallSphere();
> - /******* API ********/
> - virtual Real displacementN(){ return
> (se32.position-contactPoint).norm()-refLength;}
> - virtual Vector3r displacementT(){ relocateContactPoints();
> return contPtInTgPlane2()-contPtInTgPlane1(); }
> - virtual Real slipToDisplacementTMax(Real displacementTMax);
> - /***** end API ******/
> -
> - void setTgPlanePts(const Vector3r&, const Vector3r&);
> - void relocateContactPoints(){
> relocateContactPoints(contPtInTgPlane1(),contPtInTgPlane2()); }
> - void relocateContactPoints(const Vector3r& p1, const
> Vector3r& p2);
> -
> -
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(Dem3DofGeom_WallSphere,Dem3DofGeom,"Representation
> of contact between wall and sphere, based on Dem3DofGeom.",
> - ((Vector3r,cp1pt,,,"initial contact point on the wall,
> relative to the current contact point"))
> - ((Quaternionr,cp2rel,,,"orientation between +x and the
> reference contact point (on the sphere) in sphere-local coords"))
> - ((Real,effR2,,,"effective radius of sphere")),
> - /*ctor*/ createIndex();
> - );
> - REGISTER_CLASS_INDEX(Dem3DofGeom_WallSphere,Dem3DofGeom);
> - DECLARE_LOGGER;
> - friend class Gl1_Dem3DofGeom_WallSphere;
> - friend class Ig2_Wall_Sphere_Dem3DofGeom;
> -};
> -REGISTER_SERIALIZABLE(Dem3DofGeom_WallSphere);
> -
> -#ifdef YADE_OPENGL
> - #include<yade/pkg/common/GLDrawFunctors.hpp>
> - class Gl1_Dem3DofGeom_WallSphere:public GlIGeomFunctor{
> - public:
> - virtual void go(const shared_ptr<IGeom>&,const
> shared_ptr<Interaction>&,const shared_ptr<Body>&,const
> shared_ptr<Body>&,bool wireFrame);
> - RENDERS(Dem3DofGeom_WallSphere);
> -
> YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_Dem3DofGeom_WallSphere,GlIGeomFunctor,"Render
> interaction of wall and sphere (represented by Dem3DofGeom_WallSphere)",
> - ((bool,normal,false,,"Render interaction normal"))
> - ((bool,rolledPoints,false,,"Render points rolled
> on the spheres (tracks the original contact point)"))
> - ((bool,unrolledPoints,false,,"Render original
> contact points unrolled to the contact plane"))
> - ((bool,shear,false,,"Render shear line in the
> contact plane"))
> - ((bool,shearLabel,false,,"Render shear magnitude
> as number"))
> - );
> - };
> - REGISTER_SERIALIZABLE(Gl1_Dem3DofGeom_WallSphere);
> -#endif
> -
> -#include<yade/pkg/common/Dispatching.hpp>
> -class Ig2_Wall_Sphere_Dem3DofGeom:public IGeomFunctor{
> - public:
> - virtual bool go(const shared_ptr<Shape>& cm1, const
> shared_ptr<Shape>& cm2, const State& state1, const State& state2, const
> Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
> - virtual bool goReverse( const shared_ptr<Shape>& cm1,
> const shared_ptr<Shape>& cm2, const State& state1, const State& state2,
> const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>&
> c){
> - c->swapOrder(); return
> go(cm2,cm1,state2,state1,-shift2,force,c);
> - LOG_ERROR("!! goReverse might not work in
> Ig2_Wall_Sphere_Dem3DofGeom. IGeomDispatcher should swap interaction
> members first and call go(...) afterwards.");
> - }
> -
> - FUNCTOR2D(Wall,Sphere);
> - DEFINE_FUNCTOR_ORDER_2D(Wall,Sphere);
> -
> YADE_CLASS_BASE_DOC(Ig2_Wall_Sphere_Dem3DofGeom,IGeomFunctor,"Create/update
> contact of :yref:`Wall` and :yref:`Sphere` (:yref:`Dem3DofGeom_WallSphere`
> instance)");
> - DECLARE_LOGGER;
> -};
> -REGISTER_SERIALIZABLE(Ig2_Wall_Sphere_Dem3DofGeom);
> -
>
> === modified file 'pkg/dem/DemXDofGeom.cpp'
> --- pkg/dem/DemXDofGeom.cpp 2010-10-13 16:23:08 +0000
> +++ pkg/dem/DemXDofGeom.cpp 2013-08-29 10:30:31 +0000
> @@ -1,5 +1,3 @@
> #include"DemXDofGeom.hpp"
> -YADE_PLUGIN((GenericSpheresContact)(Dem3DofGeom));
> -Real Dem3DofGeom::displacementN(){throw;}
> -Dem3DofGeom::~Dem3DofGeom(){}
> +YADE_PLUGIN((GenericSpheresContact));
> GenericSpheresContact::~GenericSpheresContact(){}
>
> === modified file 'pkg/dem/DemXDofGeom.hpp'
> --- pkg/dem/DemXDofGeom.hpp 2011-01-24 14:45:35 +0000
> +++ pkg/dem/DemXDofGeom.hpp 2013-08-29 10:30:31 +0000
> @@ -2,12 +2,12 @@
> #pragma once
> #include<yade/core/IGeom.hpp>
>
> -/*! Abstract class that unites ScGeom and Dem3DofGeom,
> +/*! Abstract class that unites ScGeom and L3Geom,
> created for the purposes of GlobalStiffnessTimeStepper.
> It might be removed in the future. */
> class GenericSpheresContact: public IGeom{
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(GenericSpheresContact,IGeom,
> - "Class uniting :yref:`ScGeom` and :yref:`Dem3DofGeom`, for
> the purposes of :yref:`GlobalStiffnessTimeStepper`. (It might be removed
> inthe future). Do not use this class directly.",
> + "Class uniting :yref:`ScGeom` and :yref:`L3Geom`, for the
> purposes of :yref:`GlobalStiffnessTimeStepper`. (It might be removed inthe
> future). Do not use this class directly.",
> ((Vector3r,normal,,,"Unit vector oriented along the
> interaction, from particle #1, towards particle #2. |yupdate|"))
> ((Vector3r,contactPoint,,,"some reference point for the
> interaction (usually in the middle). |ycomp|"))
> ((Real,refR1,,,"Reference radius of particle #1. |ycomp|"))
> @@ -19,45 +19,3 @@
> virtual ~GenericSpheresContact(); // vtable
> };
> REGISTER_SERIALIZABLE(GenericSpheresContact);
> -
> -/*! Abstract base class for representing contact geometry of 2 elements
> that has 3 degrees of freedom:
> - * normal (1 component) and shear (Vector3r, but in plane perpendicular
> to the normal)
> - */
> -class Dem3DofGeom: public GenericSpheresContact{
> - public:
> - virtual ~Dem3DofGeom();
> -
> - // API that needs to be implemented in derived classes
> - virtual Real displacementN();
> - virtual Vector3r displacementT(){throw;}
> - virtual Real slipToDisplacementTMax(Real
> displacementTMax){throw;}; // plasticity
> - // reference radii, for contact stiffness computation; set
> to negative for nonsense values
> - // end API
> -
> - Real strainN(){
> - //if(!logCompression)
> - return displacementN()/refLength;
> - //else{Real dn=displacementN(); }
> - }
> - Vector3r strainT(){return displacementT()/refLength;}
> - Real slipToStrainTMax(Real strainTMax){return
> slipToDisplacementTMax(strainTMax*refLength)/refLength;}
> -
> -
> YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Dem3DofGeom,GenericSpheresContact,
> - "Abstract base class for representing contact
> geometry of 2 elements that has 3 degrees of freedom: normal (1 component)
> and shear (Vector3r, but in plane perpendicular to the normal).",
> - ((Real,refLength,,,"some length used to convert
> displacements to strains. |ycomp|"))
> - ((bool,logCompression,false,,"make strain go to -∞
> for length going to zero (false by default)."))
> - ((Se3r,se31,,,"Copy of body #1 se3 (needed to
> compute torque from the contact, strains etc). |yupdate|"))
> - ((Se3r,se32,,,"Copy of body #2 se3. |yupdate|")),
> - createIndex();,
> - .def("displacementN",&Dem3DofGeom::displacementN)
> - .def("displacementT",&Dem3DofGeom::displacementT)
> - .def("strainN",&Dem3DofGeom::strainN)
> - .def("strainT",&Dem3DofGeom::strainT)
> -
> .def("slipToDisplacementTMax",&Dem3DofGeom::slipToDisplacementTMax)
> -
> .def("slipToStrainTMax",&Dem3DofGeom::slipToStrainTMax)
> - );
> - REGISTER_CLASS_INDEX(Dem3DofGeom,IGeom);
> -};
> -REGISTER_SERIALIZABLE(Dem3DofGeom);
> -
> -
>
> === modified file 'pkg/dem/DomainLimiter.cpp'
> --- pkg/dem/DomainLimiter.cpp 2013-08-23 15:21:20 +0000
> +++ pkg/dem/DomainLimiter.cpp 2013-08-29 10:30:31 +0000
> @@ -81,15 +81,14 @@
> // test object types
> GenericSpheresContact*
> gsc=dynamic_cast<GenericSpheresContact*>(I->geom.get());
> ScGeom* scGeom=dynamic_cast<ScGeom*>(I->geom.get());
> - Dem3DofGeom* d3dGeom=dynamic_cast<Dem3DofGeom*>(I->geom.get());
> L3Geom* l3Geom=dynamic_cast<L3Geom*>(I->geom.get());
> L6Geom* l6Geom=dynamic_cast<L6Geom*>(I->geom.get());
> ScGeom6D* scGeom6d=dynamic_cast<ScGeom6D*>(I->geom.get());
> bool hasRot=(l6Geom || scGeom6d);
> //NormShearPhys* phys=dynamic_cast<NormShearPhys*>(I->phys.get());
> //Disabled because of warning
> if(!gsc) throw std::invalid_argument("LawTester: IGeom of
> "+strIds+" not a GenericSpheresContact.");
> - if(!scGeom && !d3dGeom && !l3Geom) throw
> std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom,
> nor Dem3DofGeom, nor L3Geom (or L6Geom).");
> - assert(!((bool)scGeom && (bool)d3dGeom && (bool)l3Geom)); //
> nonsense
> + if(!scGeom && !l3Geom) throw std::invalid_argument("LawTester:
> IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom, nor L3Geom (or
> L6Geom).");
> + assert(!((bool)scGeom && (bool)l3Geom)); // nonsense
> // get body objects
> State *state1=Body::byId(id1,scene)->state.get(),
> *state2=Body::byId(id2,scene)->state.get();
> scene->forces.sync();
> @@ -125,9 +124,7 @@
> if(scGeom6d)
> uGeom.tail<3>()=-1.*Vector3r(scGeom6d->getTwist(),scGeom6d->getBending().dot(axY),scGeom6d->getBending().dot(axZ));
> }
> else{ // d3dGeom
> - throw runtime_error("LawTester:
> Dem3DofGeom not yet supported.");
> - // essentially copies code from ScGeom,
> which is not very nice indeed; oh well…
> - // Vector3r
> vRel=(state2->vel+state2->angVel.cross(-gsc->refR2*gsc->normal))-(state1->vel+state1->angVel.cross(gsc->refR1*gsc->normal));
> + throw runtime_error("Geom type not yet
> supported.");
> }
> }
> // update the transformation
>
> === modified file 'pkg/dem/FlowEngine.cpp'
> --- pkg/dem/FlowEngine.cpp 2013-07-26 18:16:04 +0000
> +++ pkg/dem/FlowEngine.cpp 2013-08-30 13:27:06 +0000
> @@ -21,6 +21,10 @@
> #include <boost/date_time.hpp>
> #include <boost/date_time/posix_time/posix_time.hpp>
>
> +#ifdef LINSOLV
> +#include <cholmod.h>
> +#endif
> +
> #include "FlowEngine.hpp"
>
> CREATE_LOGGER ( FlowEngine );
> @@ -219,6 +223,20 @@
> flow->x_min = 1000.0, flow->x_max = -10000.0, flow->y_min =
> 1000.0, flow->y_max = -10000.0, flow->z_min = 1000.0, flow->z_max =
> -10000.0;
> }
>
> +#ifdef LINSOLV
> +template<class Solver>
> +void FlowEngine::setForceMetis ( Solver& flow, bool force )
> +{
> + if (force) {
> + flow->eSolver.cholmod().nmethods=1;
> + flow->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;
> + } else cholmod_defaults(&(flow->eSolver.cholmod()));
> +}
> +
> +template<class Solver>
> +bool FlowEngine::getForceMetis ( Solver& flow ) {return
> (flow->eSolver.cholmod().nmethods==1);}
> +#endif
> +
> template<class Solver>
> void FlowEngine::Build_Triangulation ( Solver& flow )
> {
>
> === modified file 'pkg/dem/FlowEngine.hpp'
> --- pkg/dem/FlowEngine.hpp 2013-08-22 18:35:55 +0000
> +++ pkg/dem/FlowEngine.hpp 2013-08-30 13:07:40 +0000
> @@ -59,6 +59,10 @@
> int ReTrg;
> int ellapsedIter;
> TPL void initSolver (Solver& flow);
> + #ifdef LINSOLV
> + TPL void setForceMetis (Solver& flow, bool force);
> + TPL bool getForceMetis (Solver& flow);
> + #endif
> TPL void Triangulate (Solver& flow);
> TPL void AddBoundary (Solver& flow);
> TPL void Build_Triangulation (double P_zero, Solver& flow);
> @@ -155,6 +159,13 @@
> #ifdef LINSOLV
> void _exportMatrix(string filename)
> {exportMatrix(filename,solver);}
> void _exportTriplets(string filename)
> {exportTriplets(filename,solver);}
> + void _setForceMetis (bool force)
> {setForceMetis(solver,force);}
> + bool _getForceMetis () {return getForceMetis
> (solver);}
> + void cholmodStats() {
> + cerr << cholmod_print_common("PFV
> Cholmod factorization",&(solver->eSolver.cholmod()))<<endl;
> + cerr << "cholmod method:" <<
> solver->eSolver.cholmod().selected<<endl;
> + cerr << "METIS
> called:"<<solver->eSolver.cholmod().called_nd<<endl;}
> + bool metisUsed() {return
> bool(solver->eSolver.cholmod().called_nd);}
> #endif
> python::list _getConstrictions(bool all) {return
> getConstrictions(all,solver);}
> python::list _getConstrictionsFull(bool all) {return
> getConstrictionsFull(all,solver);}
> @@ -164,7 +175,7 @@
> virtual void action();
> virtual void backgroundAction();
>
> -
> YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(FlowEngine,PartialEngine,"An
> engine to solve flow problem in saturated granular media. Model description
> can be found in [Chareyre2012a]_ and [Catalano2013a]_. See the example
> script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n..
> note::\n\t Multi-threading seems to work fine for Cholesky decomposition,
> but it fails for the solve phase in which -j1 is the fastest, here we
> specify thread numbers independently using
> :yref:`FlowEngine::numFactorizeThreads` and
> :yref:`FlowEngine::numSolveThreads`. These multhreading settings are only
> impacting the behaviour of openblas library and are relatively independant
> of :yref:`FlowEngine::multithread`. However, the settings have to be
> globally consistent. For instance, :yref:`FlowEngine::multithread`=True
> with multithread
> :yref:`FlowEngine::numSolveThreads`=:yref:`FlowEngine::numFactorizeThreads`=4
> means that openblas will mobilize 8 processors at some point, if the system
> doesn't have so many procs. it will hurt performance.",
> +
> YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(FlowEngine,PartialEngine,"An
> engine to solve flow problem in saturated granular media. Model description
> can be found in [Chareyre2012a]_ and [Catalano2013a]_. See the example
> script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n..
> note::Multi-threading seems to work fine for Cholesky decomposition, but it
> fails for the solve phase in which -j1 is the fastest, here we specify
> thread numbers independently using :yref:`FlowEngine::numFactorizeThreads`
> and :yref:`FlowEngine::numSolveThreads`. These multhreading settings are
> only impacting the behaviour of openblas library and are relatively
> independant of :yref:`FlowEngine::multithread`. However, the settings have
> to be globally consistent. For instance,
> :yref:`multithread<FlowEngine::multithread>` =True with
> yref:`numFactorizeThreads<FlowEngine::numFactorizeThreads>` =
> yref:`numSolveThreads<FlowEngine::numSolveThreads>` = 4 implies that
> openblas will mobilize 8 processors at some point. If the system does not
> have so many procs. it will hurt performance.",
>
> ((bool,isActivated,true,,"Activates Flow Engine"))
> ((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."))
> @@ -188,7 +199,7 @@
> ((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,permeabilityFactor,0.0,,"permability multiplier"))
> +
> ((double,permeabilityFactor,1.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"))
> @@ -217,6 +228,7 @@
> #ifdef EIGENSPARSE_LIB
> ((int, numSolveThreads, 1,,"number
> of openblas threads in the solve phase."))
> ((int, numFactorizeThreads,
> 1,,"number of openblas threads in the factorization phase"))
> +// ((bool, forceMetis, 0,,"If true,
> METIS is used for matrix preconditioning, else Cholmod is free to choose
> the best method (which may be METIS to, depending on the matrix). See
> ``nmethods`` in Cholmod documentation"))
> #endif
> ,
> /*deprec*/
> @@ -240,7 +252,7 @@
>
> .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.\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("getBoundaryFlux",&FlowEngine::_getBoundaryFlux,(python::arg("boundary")),"Get
> total flux through boundary defined by its body id.\n\n.. note:: 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.")
> @@ -261,6 +273,9 @@
> #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.")
> +
> .def("cholmodStats",&FlowEngine::cholmodStats,"get statistics of cholmod
> solver activity")
> +
> .def("metisUsed",&FlowEngine::metisUsed,"check wether metis lib is
> effectively used")
> +
> .add_property("forceMetis",&FlowEngine::_getForceMetis,&FlowEngine::_setForceMetis,"If
> true, METIS is used for matrix preconditioning, else Cholmod is free to
> choose the best method (which may be METIS to, depending on the matrix).
> See ``nmethods`` in Cholmod documentation")
> #endif
> )
> DECLARE_LOGGER;
>
> === modified file 'pkg/dem/HertzMindlin.cpp'
> --- pkg/dem/HertzMindlin.cpp 2012-02-16 16:05:15 +0000
> +++ pkg/dem/HertzMindlin.cpp 2013-09-07 19:40:12 +0000
> @@ -31,8 +31,8 @@
>
> void Ip2_FrictMat_FrictMat_MindlinPhys::go(const shared_ptr<Material>&
> b1,const shared_ptr<Material>& b2, const shared_ptr<Interaction>&
> interaction){
> if(interaction->phys) return; // no updates of an already existing
> contact necessary
> - shared_ptr<MindlinPhys> mindlinPhys(new MindlinPhys());
> - interaction->phys = mindlinPhys;
> + shared_ptr<MindlinPhys> contactPhysics(new MindlinPhys());
> + interaction->phys = contactPhysics;
> FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
> FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
>
> @@ -64,20 +64,20 @@
> Real Rmean = (Da+Db)/2.; // mean radius
> Real Kno = 4./3.*E*sqrt(R); // coefficient for normal stiffness
> Real Kso = 2*sqrt(4*R)*G/(2-V); // coefficient for shear stiffness
> - Real frictionAngle = std::min(fa,fb);
> + Real frictionAngle = (!frictAngle) ? std::min(fa,fb) :
> (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle);
>
> Real Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force
> as predicted by DMT theory
>
> /* pass values calculated from above to MindlinPhys */
> - mindlinPhys->tangensOfFrictionAngle = std::tan(frictionAngle);
> - //mindlinPhys->prevNormal = scg->normal; // used to compute
> relative rotation
> - mindlinPhys->kno = Kno; // this is just a coeff
> - mindlinPhys->kso = Kso; // this is just a coeff
> - mindlinPhys->adhesionForce = Adhesion;
> + contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
> + //contactPhysics->prevNormal = scg->normal; // used to compute
> relative rotation
> + contactPhysics->kno = Kno; // this is just a coeff
> + contactPhysics->kso = Kso; // this is just a coeff
> + contactPhysics->adhesionForce = Adhesion;
>
> - mindlinPhys->kr = krot;
> - mindlinPhys->ktw = ktwist;
> - mindlinPhys->maxBendPl = eta*Rmean; // does this make sense? why
> do we take Rmean?
> + contactPhysics->kr = krot;
> + contactPhysics->ktw = ktwist;
> + contactPhysics->maxBendPl = eta*Rmean; // does this make sense?
> why do we take Rmean?
>
> /* compute viscous coefficients */
> if(en && betan) throw
> std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinPhys: only one of en,
> betan can be specified.");
> @@ -86,13 +86,13 @@
> // en or es specified, just compute alpha, otherwise alpha remains
> 0
> if(en || es){
> Real logE = log((*en)(mat1->id,mat2->id));
> - mindlinPhys->alpha =
> -sqrt(5/6.)*2*logE/sqrt(pow(logE,2)+pow(Mathr::PI,2))*sqrt(2*E*sqrt(R)); //
> (see Tsuji, 1992)
> + contactPhysics->alpha =
> -sqrt(5/6.)*2*logE/sqrt(pow(logE,2)+pow(Mathr::PI,2))*sqrt(2*E*sqrt(R)); //
> (see Tsuji, 1992)
> }
>
> // betan specified, use that value directly; otherwise give zero
> else{
> - mindlinPhys->betan=betan ? (*betan)(mat1->id,mat2->id) : 0;
> - mindlinPhys->betas=betas ? (*betas)(mat1->id,mat2->id) :
> mindlinPhys->betan;
> + contactPhysics->betan=betan ? (*betan)(mat1->id,mat2->id)
> : 0;
> + contactPhysics->betas=betas ? (*betas)(mat1->id,mat2->id)
> : contactPhysics->betan;
> }
> }
>
>
> === modified file 'pkg/dem/HertzMindlin.hpp'
> --- pkg/dem/HertzMindlin.hpp 2013-05-29 11:54:22 +0000
> +++ pkg/dem/HertzMindlin.hpp 2013-09-07 19:40:12 +0000
> @@ -86,6 +86,7 @@
> ((shared_ptr<MatchMaker>,es,,,"Shear coefficient
> of restitution $e_s$."))
> ((shared_ptr<MatchMaker>,betan,,,"Normal viscous
> damping coefficient $\\beta_n$."))
> ((shared_ptr<MatchMaker>,betas,,,"Shear viscous
> damping coefficient $\\beta_s$."))
> + ((shared_ptr<MatchMaker>,frictAngle,,,"Instance of
> :yref:`MatchMaker` determining how to compute the friction angle of an
> interaction. If ``None``, minimum value is used."))
> );
> DECLARE_LOGGER;
> };
>
> === modified file 'pkg/dem/Shop.cpp'
> --- pkg/dem/Shop.cpp 2013-09-08 11:12:42 +0000
> +++ pkg/dem/Shop.cpp 2013-09-08 11:19:06 +0000
> @@ -581,44 +581,23 @@
> Vector3r normalStress,shearStress;
> if(!I->isReal()) continue;
>
> - //NormShearPhys + Dem3DofGeom
> - const NormShearPhys* physNSP =
> YADE_CAST<NormShearPhys*>(I->phys.get());
> - //Dem3DofGeom*
> geom=YADE_CAST<Dem3DofGeom*>(I->geom.get()); //For the moment only for
> Dem3DofGeom!!!
> - // FIXME: slower, but does not crash
> - Dem3DofGeom*
> geomDDG=dynamic_cast<Dem3DofGeom*>(I->geom.get()); //For the moment only
> for Dem3DofGeom!!!
> -
> - //FrictPhys + ScGeom
> const FrictPhys* physFP =
> YADE_CAST<FrictPhys*>(I->phys.get());
> - ScGeom* geomScG=dynamic_cast<ScGeom*>(I->geom.get());
> + ScGeom* geomScG=YADE_CAST<ScGeom*>(I->geom.get());
>
> const Body::id_t id1=I->getId1(), id2=I->getId2();
>
> - if(((physNSP) and (geomDDG)) or ((physFP) and (geomScG))){
> - if ((physNSP) and (geomDDG)) {
> - Real
> minRad=(geomDDG->refR1<=0?geomDDG->refR2:(geomDDG->refR2<=0?geomDDG->refR1:min(geomDDG->refR1,geomDDG->refR2)));
> - Real crossSection=Mathr::PI*pow(minRad,2);
> -
> -
> normalStress=((1./crossSection)*geomDDG->normal.dot(physNSP->normalForce))*geomDDG->normal;
> - for(int i=0; i<3; i++){
> - int ix1=(i+1)%3,ix2=(i+2)%3;
> -
> shearStress[i]=geomDDG->normal[ix1]*physNSP->shearForce[ix1]+geomDDG->normal[ix2]*physNSP->shearForce[ix2];
> - shearStress[i]/=crossSection;
> - }
> - } else if ((physFP) and (geomScG)) {
> - Real
> minRad=(geomScG->radius1<=0?geomScG->radius2:(geomScG->radius2<=0?geomScG->radius1:min(geomScG->radius1,geomScG->radius2)));
> - Real crossSection=Mathr::PI*pow(minRad,2);
> + if((physFP) and (geomScG)){
> + Real
> minRad=(geomScG->radius1<=0?geomScG->radius2:(geomScG->radius2<=0?geomScG->radius1:min(geomScG->radius1,geomScG->radius2)));
> + Real crossSection=Mathr::PI*pow(minRad,2);
>
> -
> normalStress=((1./crossSection)*geomScG->normal.dot(physFP->normalForce))*geomScG->normal;
> - for(int i=0; i<3; i++){
> - int ix1=(i+1)%3,ix2=(i+2)%3;
> -
> shearStress[i]=geomScG->normal[ix1]*physFP->shearForce[ix1]+geomScG->normal[ix2]*physFP->shearForce[ix2];
> - shearStress[i]/=crossSection;
> - }
> +
> normalStress=((1./crossSection)*geomScG->normal.dot(physFP->normalForce))*geomScG->normal;
> + for(int i=0; i<3; i++){
> + int ix1=(i+1)%3,ix2=(i+2)%3;
> +
> shearStress[i]=geomScG->normal[ix1]*physFP->shearForce[ix1]+geomScG->normal[ix2]*physFP->shearForce[ix2];
> + shearStress[i]/=crossSection;
> }
> -
> bodyStates[id1].normStress+=normalStress;
> bodyStates[id2].normStress+=normalStress;
> -
> bodyStates[id1].shearStress+=shearStress;
> bodyStates[id2].shearStress+=shearStress;
> }
> @@ -818,33 +797,6 @@
> return stressTensor/volume;
> }
>
> -/*
> -Matrix3r Shop::stressTensorOfPeriodicCell(bool smallStrains){
> - Scene* scene=Omega::instance().getScene().get();
> - if (!scene->isPeriodic){ throw runtime_error("Can't compute stress
> of periodic cell in aperiodic simulation."); }
> -// if (smallStrains){volume =
> scene->getVolume()cell->refSize[0]*scene->cell->refSize[1]*scene->cell->refSize[2];}
> -// else volume =
> scene->cell->trsf.determinant()*scene->cell->refSize[0]*scene->cell->refSize[1]*scene->cell->refSize[2];
> - // Using the function provided by Cell (BC)
> - Real volume = scene->cell->getVolume();
> - Matrix3r stress = Matrix3r::Zero();
> - FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
> - if(!I->isReal()) continue;
> - Dem3DofGeom* geom=YADE_CAST<Dem3DofGeom*>(I->geom.get());
> - NormShearPhys*
> phys=YADE_CAST<NormShearPhys*>(I->phys.get());
> - Real l;
> - if (smallStrains){l = geom->refLength;}
> - else l=(geom->se31.position-geom->se32.position).norm();
> - Vector3r& n=geom->normal;
> - Vector3r& fT=phys->shearForce;
> - Real fN=phys->normalForce.dot(n);
> -
> - stress += l*(fN*n*n.transpose() + .5*(fT*n.transpose() +
> n*fT.transpose()));
> - }
> - stress/=volume;
> - return stress;
> -}
> -*/
> -
> void Shop::getStressLWForEachBody(vector<Matrix3r>& bStresses){
> const shared_ptr<Scene>& scene=Omega::instance().getScene();
> bStresses.resize(scene->bodies->size());
>
> === modified file 'pkg/dem/Shop.hpp'
> --- pkg/dem/Shop.hpp 2013-09-08 11:12:42 +0000
> +++ pkg/dem/Shop.hpp 2013-09-08 11:19:06 +0000
> @@ -127,8 +127,6 @@
> static Matrix3r getStress(Real volume=0);
> static Matrix3r getCapillaryStress(Real volume=0);
> static Matrix3r stressTensorOfPeriodicCell() {
> LOG_WARN("Shop::stressTensorOfPeriodicCelli is DEPRECATED: use getStress
> instead"); return Shop::getStress(); }
> - //! This version is restricted to periodic BCs and Dem3Dof
> - static Matrix3r stressTensorOfPeriodicCell(bool
> smallStrains=true);
> //! Compute overall ("macroscopic") stress of periodic
> cell, returning 2 tensors
> //! (contribution of normal and shear forces)
> static py::tuple normalShearStressTensors(bool
> compressionPositive=false, bool splitNormalTensor=false, Real
> thresholdForce=NaN);
>
> === modified file 'pkg/dem/VTKRecorder.hpp'
> --- pkg/dem/VTKRecorder.hpp 2013-04-24 19:49:00 +0000
> +++ pkg/dem/VTKRecorder.hpp 2013-08-29 10:30:31 +0000
> @@ -23,7 +23,7 @@
> ((bool,multiblock,false,,"Use multi-block
> (``.vtm``) files to store data, rather than separate ``.vtu`` files."))
> #endif
> ((string,fileName,"",,"Base file name; it will be appended
> with {spheres,intrs,facets}-243100.vtu (unless *multiblock* is ``True``)
> depending on active recorders and step number (243100 in this case). It can
> contain slashes, but the directory must exist already."))
> -
> ((vector<string>,recorders,vector<string>(1,string("all")),,"List of active
> recorders (as strings). ``all`` (the default value) enables all base and
> generic recorders.\n\n.. admonition:: Base recorders\n\n\tBase recorders
> save the geometry (unstructured grids) on which other data is defined. They
> are implicitly activated by many of the other recorders. Each of them
> creates a new file (or a block, if :yref:`multiblock
> <VTKRecorder.multiblock>` is set).\n\n\t``spheres``\n\t\tSaves positions
> and radii (``radii``) of :yref:`spherical<Sphere>`
> particles.\n\t``facets``\n\t\tSave :yref:`facets<Facet>` positions
> (vertices).\n\t``boxes``\n\t\tSave :yref:`boxes<Box>` positions
> (edges).\n\t``intr``\n\t\tStore interactions as lines between nodes at
> respective particles positions. Additionally stores magnitude of normal
> (``forceN``) and shear (``absForceT``) forces on interactions (the
> :yref:`geom<Interaction.geom> must be of type :yref:`NormShearPhys`).
> \n\n.. admonition:: Generic recorders\n\n\tGeneric recorders do not depend
> on specific model being used and save commonly useful
> data.\n\n\t``id``\n\t\tSaves id's (field ``id``) of spheres; active only if
> ``spheres`` is active.\n\t``mass``\n\t\tSaves masses (field ``mass``) of
> spheres; active only if ``spheres`` is active.\n\t``clumpId``\n\t\tSaves
> id's of clumps to which each sphere belongs (field ``clumpId``); active
> only if ``spheres`` is active.\n\t``colors``\n\t\tSaves colors of
> :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``);
> only active if ``spheres`` or ``facets`` are
> activated.\n\t``mask``\n\t\tSaves groupMasks of :yref:`spheres<Sphere>` and
> of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or
> ``facets`` are activated.\n\t``materialId``\n\t\tSaves materialID of
> :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if
> ``spheres`` or ``facets`` are activated.\n\t``velocity``\n\t\tSaves linear
> and angular velocities of spherical particles as Vector3 and length(fields
> ``linVelVec``, ``l
> ...
>
> [Zpráva zkrácena]
References
