yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11449
[Branch ~yade-pkg/yade/git-trunk] Rev 3444: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Jérôme Duriez (jduriez)
------------------------------------------------------------
revno: 3444 [merge]
committer: cyuan <chaoyuan2012@xxxxxxxxx>
timestamp: Thu 2014-07-17 12:21:07 +0800
message:
Merge github.com:yade/trunk into chaoUnsat
removed:
pkg/dem/CohFrictMat.cpp
pkg/dem/CohFrictMat.hpp
pkg/dem/CohFrictPhys.cpp
pkg/dem/CohFrictPhys.hpp
pkg/dem/InelastCohFrictMat.cpp
pkg/dem/InelastCohFrictMat.hpp
pkg/dem/InelastCohFrictPhys.cpp
pkg/dem/InelastCohFrictPhys.hpp
pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp
pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp
pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp
pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp
pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp
pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp
pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp
pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.hpp
pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.cpp
pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp
pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp
pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp
pkg/dem/KinemCNDEngine.cpp
pkg/dem/KinemCNDEngine.hpp
pkg/dem/KinemCNLEngine.cpp
pkg/dem/KinemCNLEngine.hpp
pkg/dem/KinemCNSEngine.cpp
pkg/dem/KinemCNSEngine.hpp
pkg/dem/KinemCTDEngine.cpp
pkg/dem/KinemCTDEngine.hpp
pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp
pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp
pkg/dem/NormalInelasticMat.cpp
pkg/dem/NormalInelasticMat.hpp
pkg/dem/NormalInelasticityLaw.cpp
pkg/dem/NormalInelasticityLaw.hpp
pkg/dem/NormalInelasticityPhys.cpp
pkg/dem/NormalInelasticityPhys.hpp
added:
pkg/dem/InelastCohFrictPM.cpp
pkg/dem/InelastCohFrictPM.hpp
pkg/dem/KinemC__Engine.cpp
pkg/dem/KinemC__Engine.hpp
pkg/dem/NormalInelasticPM.cpp
pkg/dem/NormalInelasticPM.hpp
modified:
CMakeLists.txt
LICENSE
core/Omega.cpp
core/Omega.hpp
core/main/pyboot.cpp
doc/sphinx/user.rst
gui/qt4/GLViewer.cpp
gui/qt4/GLViewer.hpp
gui/qt4/GLViewerDisplay.cpp
gui/qt4/GLViewerMouse.cpp
lib/base/Math.hpp
lib/factory/DynLibManager.cpp
lib/serialization/Serializable.hpp
lib/triangulation/FlowBoundingSphere.ipp
lib/triangulation/Network.ipp
pkg/common/Cylinder.hpp
pkg/common/Grid.hpp
pkg/common/PersistentTriangulationCollider.hpp
pkg/dem/CapillaryPhys.cpp
pkg/dem/CapillaryPhys.hpp
pkg/dem/CapillaryTriaxialTest.cpp
pkg/dem/CohesiveFrictionalContactLaw.cpp
pkg/dem/CohesiveFrictionalContactLaw.hpp
pkg/dem/CohesiveTriaxialTest.cpp
pkg/dem/ConcretePM.hpp
pkg/dem/ElasticContactLaw.hpp
pkg/dem/FrictPhys.cpp
pkg/dem/FrictPhys.hpp
pkg/dem/FrictViscoPM.hpp
pkg/dem/HertzMindlin.cpp
pkg/dem/HertzMindlin.hpp
pkg/dem/Ip2_ElastMat.cpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp
pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp
pkg/dem/MicroMacroAnalyser.hpp
pkg/dem/Shop_01.cpp
pkg/dem/Shop_02.cpp
pkg/dem/SimpleShear.cpp
pkg/dem/TriaxialTest.cpp
pkg/dem/VTKRecorder.cpp
pkg/dem/VTKRecorder.hpp
pkg/lbm/HydrodynamicsLawLBM.cpp
py/mathWrap/miniEigen.cpp
py/wrapper/yadeWrapper.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'CMakeLists.txt'
--- CMakeLists.txt 2014-07-09 07:07:14 +0000
+++ CMakeLists.txt 2014-07-17 04:21:07 +0000
@@ -14,8 +14,8 @@
# ENABLE_GL2PS: enable GL2PS-option (ON by default)
# ENABLE_LINSOLV: enable LINSOLV-option (ON by default)
# ENABLE_PFVFLOW: enable PFVFLOW-option, FlowEngine (ON by default)
+# ENABLE_LBMFLOW: enable LBMFLOW-option, LBM_ENGINE (ON by default)
# ENABLE_SPH: enable SPH-option, Smoothed Particle Hydrodynamics (OFF by default, experimental)
-# ENABLE_LBMFLOW: enable LBMFLOW-option, LBM_ENGINE (OFF by default)
# ENABLE_LIQMIGRATION: enable LIQMIGRATION-option, see [Mani2013] for details (OFF by default)
# ENABLE_MASK_ARBITRARY: enable MASK_ARBITRARY option (OFF by default)
# runtimePREFIX: used for packaging, when install directory is not the same is runtime directory (/usr/local by default)
@@ -63,13 +63,7 @@
#===========================================================
-ADD_DEFINITIONS(" -DYADE_PTR_CAST=boost::static_pointer_cast -DYADE_CAST=static_cast -DYADE_PTR_DYN_CAST=boost::dynamic_pointer_cast ")
-IF (CXX11)
- SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x")
- ADD_DEFINITIONS(" -DCXX11 -DTYPEOF=decltype")
-ELSE (CXX11)
- ADD_DEFINITIONS(" -DTYPEOF=typeof")
-ENDIF (CXX11)
+SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x")
#===========================================================
@@ -134,7 +128,7 @@
OPTION(ENABLE_LINSOLV "Enable direct solver for the flow engines (experimental)" ${DEFAULT})
OPTION(ENABLE_PFVFLOW "Enable flow engine (experimental)" ${DEFAULT})
OPTION(ENABLE_SPH "Enable SPH" OFF)
-OPTION(ENABLE_LBMFLOW "Enable LBM engine (very experimental)" OFF)
+OPTION(ENABLE_LBMFLOW "Enable LBM engine (very experimental)" ON)
OPTION(ENABLE_LIQMIGRATION "Enable liquid control (very experimental), see [Mani2013] for details" OFF)
OPTION(ENABLE_MASK_ARBITRARY "Enable arbitrary precision of bitmask variables (only Body::groupMask yet implemented) (experimental). Use -DMASK_ARBITRARY_SIZE=int to set number of used bits (256 by default)" OFF)
@@ -149,7 +143,6 @@
MESSAGE(FATAL_ERROR "Minimal required Eigen3 version is 3.2.1, please update Eigen3!")
ENDIF ((${EIGEN3_MAJOR_VERSION} LESS 2) OR ((${EIGEN3_MAJOR_VERSION} EQUAL 2) AND (${EIGEN3_MINOR_VERSION} LESS 1)))
- SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} Eigen3")
IF (NOT VECTORIZE)
MESSAGE(STATUS "Disable vectorization")
ADD_DEFINITIONS("-DEIGEN_DONT_VECTORIZE -DEIGEN_DONT_ALIGN -DEIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT")
@@ -274,9 +267,9 @@
IF(ENABLE_PFVFLOW)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DFLOW_ENGINE")
- SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} PFVflow")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} PFVFLOW")
ELSE(ENABLE_PFVFLOW)
- SET(DISABLED_FEATS "${DISABLED_FEATS} PFVflow")
+ SET(DISABLED_FEATS "${DISABLED_FEATS} PFVFLOW")
ENDIF(ENABLE_PFVFLOW)
ELSE(CGAL_FOUND AND GMP_FOUND AND (NOT("${CMAKE_CXX_COMPILER} ${CMAKE_CXX_COMPILER_ARG1}" MATCHES ".*clang")))
@@ -284,7 +277,7 @@
SET(DISABLED_FEATS "${DISABLED_FEATS} CGAL")
SET(ENABLE_CGAL OFF)
IF(ENABLE_PFVFLOW)
- SET(DISABLED_FEATS "${DISABLED_FEATS} PFVflow")
+ SET(DISABLED_FEATS "${DISABLED_FEATS} PFVFLOW")
MESSAGE(STATUS "CGAL NOT found: PFVFLOW disabled")
ENDIF(ENABLE_PFVFLOW)
@@ -304,14 +297,14 @@
MESSAGE(STATUS "Found Cholmod")
MESSAGE(STATUS "Found OpenBlas")
MESSAGE(STATUS "Found Metis")
- SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} LinSolv")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} LINSOLV")
ELSE(CHOLMOD_FOUND AND OPENBLAS_FOUND AND METIS_FOUND)
MESSAGE(STATUS "Missing dependency for LINSOLV, disabled")
- SET(DISABLED_FEATS "${DISABLED_FEATS} LinSolv")
+ SET(DISABLED_FEATS "${DISABLED_FEATS} LINSOLV")
SET(ENABLE_LINSOLV OFF)
ENDIF(CHOLMOD_FOUND AND OPENBLAS_FOUND AND METIS_FOUND)
ELSE(ENABLE_LINSOLV)
- SET(DISABLED_FEATS "${DISABLED_FEATS} LinSolv")
+ SET(DISABLED_FEATS "${DISABLED_FEATS} LINSOLV")
ENDIF(ENABLE_LINSOLV)
#===============================================
IF(ENABLE_SPH)
=== modified file 'LICENSE'
--- LICENSE 2010-01-07 12:55:44 +0000
+++ LICENSE 2014-07-14 20:08:34 +0000
@@ -1,12 +1,12 @@
- GNU GENERAL PUBLIC LICENSE
- Version 2, June 1991
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
- Copyright (C) 1989, 1991 Free Software Foundation, Inc.
- 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
- Preamble
+ Preamble
The licenses for most software are designed to take away your
freedom to share and change it. By contrast, the GNU General Public
@@ -15,7 +15,7 @@
General Public License applies to most of the Free Software
Foundation's software and to any other program whose authors commit to
using it. (Some other Free Software Foundation software is covered by
-the GNU Library General Public License instead.) You can apply it to
+the GNU Lesser General Public License instead.) You can apply it to
your programs, too.
When we speak of free software, we are referring to freedom, not
@@ -55,8 +55,8 @@
The precise terms and conditions for copying, distribution and
modification follow.
-
- GNU GENERAL PUBLIC LICENSE
+
+ GNU GENERAL PUBLIC LICENSE
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
0. This License applies to any program or other work which contains
@@ -110,7 +110,7 @@
License. (Exception: if the Program itself is interactive but
does not normally print such an announcement, your work based on
the Program is not required to print an announcement.)
-
+
These requirements apply to the modified work as a whole. If
identifiable sections of that work are not derived from the Program,
and can be reasonably considered independent and separate works in
@@ -168,7 +168,7 @@
access to copy the source code from the same place counts as
distribution of the source code, even though third parties are not
compelled to copy the source along with the object code.
-
+
4. You may not copy, modify, sublicense, or distribute the Program
except as expressly provided under this License. Any attempt
otherwise to copy, modify, sublicense or distribute the Program is
@@ -225,7 +225,7 @@
This section is intended to make thoroughly clear what is believed to
be a consequence of the rest of this License.
-
+
8. If the distribution and/or use of the Program is restricted in
certain countries either by patents or by copyrighted interfaces, the
original copyright holder who places the Program under this License
@@ -255,7 +255,7 @@
of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally.
- NO WARRANTY
+ NO WARRANTY
11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
@@ -277,9 +277,9 @@
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
- END OF TERMS AND CONDITIONS
-
- How to Apply These Terms to Your New Programs
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make it
@@ -303,17 +303,16 @@
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
- You should have received a copy of the GNU General Public License
- along with this program; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
Also add information on how to contact you by electronic and paper mail.
If the program is interactive, make it output a short notice like this
when it starts in an interactive mode:
- Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision version 69, Copyright (C) year name of author
Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type `show c' for details.
@@ -336,11 +335,5 @@
This General Public License does not permit incorporating your program into
proprietary programs. If your program is a subroutine library, you may
consider it more useful to permit linking proprietary applications with the
-library. If this is what you want to do, use the GNU Library General
+library. If this is what you want to do, use the GNU Lesser General
Public License instead of this License.
-
-13. As a special exception, you have permission to link this program
-with the CGAL library and distribute executables, as long as you
-follow the requirements of the GNU GPL in regard to all of the
-software in the executable aside from CGAL.
-
=== modified file 'core/Omega.cpp'
--- core/Omega.cpp 2014-07-04 17:46:32 +0000
+++ core/Omega.cpp 2014-07-14 20:08:34 +0000
@@ -15,9 +15,6 @@
#include<yade/lib/base/Math.hpp>
#include<yade/lib/multimethods/FunctorWrapper.hpp>
#include<yade/lib/multimethods/Indexable.hpp>
-#include<boost/filesystem/operations.hpp>
-#include<boost/filesystem/convenience.hpp>
-#include<boost/filesystem/exception.hpp>
#include<boost/algorithm/string.hpp>
#include<boost/thread/mutex.hpp>
@@ -245,7 +242,7 @@
RenderMutexLock lock;
if(isMem){
istringstream iss(memSavedSimulations[f]);
- yade::ObjectIO::load<TYPEOF(scene),boost::archive::binary_iarchive>(iss,"scene",scene);
+ yade::ObjectIO::load<decltype(scene),boost::archive::binary_iarchive>(iss,"scene",scene);
} else {
yade::ObjectIO::load(f,"scene",scene);
}
@@ -267,7 +264,7 @@
if(boost::algorithm::starts_with(f,":memory:")){
if(memSavedSimulations.count(f)>0 && !quiet) LOG_INFO("Overwriting in-memory saved simulation "<<f);
ostringstream oss;
- yade::ObjectIO::save<TYPEOF(scene),boost::archive::binary_oarchive>(oss,"scene",scene);
+ yade::ObjectIO::save<decltype(scene),boost::archive::binary_oarchive>(oss,"scene",scene);
memSavedSimulations[f]=oss.str();
}
else {
=== modified file 'core/Omega.hpp'
--- core/Omega.hpp 2014-07-04 17:46:32 +0000
+++ core/Omega.hpp 2014-07-14 20:08:34 +0000
@@ -13,16 +13,14 @@
#include <Python.h>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <fstream>
-#include <set>
-#include <list>
#include <time.h>
#include <boost/thread/thread.hpp>
#include <iostream>
-#include<yade/lib/base/Math.hpp>
-#include<yade/lib/factory/ClassFactory.hpp>
+#include <yade/lib/base/Math.hpp>
+#include <yade/lib/factory/ClassFactory.hpp>
-#include<yade/lib/base/Singleton.hpp>
+#include <yade/lib/base/Singleton.hpp>
#include "SimulationFlow.hpp"
=== modified file 'core/main/pyboot.cpp'
--- core/main/pyboot.cpp 2014-07-03 17:20:40 +0000
+++ core/main/pyboot.cpp 2014-07-14 20:08:34 +0000
@@ -3,8 +3,6 @@
#include<signal.h>
-#include<boost/filesystem/convenience.hpp>
-
#ifdef YADE_DEBUG
void crashHandler(int sig){
switch(sig){
=== modified file 'doc/sphinx/user.rst'
--- doc/sphinx/user.rst 2014-05-06 15:32:52 +0000
+++ doc/sphinx/user.rst 2014-07-10 09:57:20 +0000
@@ -495,7 +495,7 @@
Individual interactions on demand
----------------------------------
-It is possible to create an interaction between a pair of particles independently of collision detection using :yref:`createInteraction<yade._utils.createInteraction>`. This function looks for and uses matching ``Ig2`` and ``Ip2`` functors. Interaction will be created regardless of distance between given particles (by passing a special parameter to the ``Ig2`` functor to force creation of the interaction even without any geometrical contact). Appropriate constitutive law should be used to avoid deletion of the interaction at the next simulation step.
+It is possible to create an interaction between a pair of particles independently of collision detection using :yref:`createInteraction<yade.utils.createInteraction>`. This function looks for and uses matching ``Ig2`` and ``Ip2`` functors. Interaction will be created regardless of distance between given particles (by passing a special parameter to the ``Ig2`` functor to force creation of the interaction even without any geometrical contact). Appropriate constitutive law should be used to avoid deletion of the interaction at the next simulation step.
.. ipython::
@@ -523,7 +523,7 @@
# created by functors in InteractionLoop
Yade [2]: i.geom, i.phys
-This method will be rather slow if many interaction are to be created (the functor lookup will be repeated for each of them). In such case, ask on yade-dev@xxxxxxxxxxxxxxxxxxx to have the :yref:`createInteraction<yade._utils.createInteraction>` function accept list of pairs id's as well.
+This method will be rather slow if many interaction are to be created (the functor lookup will be repeated for each of them). In such case, ask on yade-dev@xxxxxxxxxxxxxxxxxxx to have the :yref:`createInteraction<yade.utils.createInteraction>` function accept list of pairs id's as well.
Base engines
=============
@@ -919,7 +919,8 @@
A special engine :yref:`PyRunner` can be used to periodically call python code, specified via the ``command`` parameter. Periodicity can be controlled by specifying computation time (``realPeriod``), virutal time (``virtPeriod``) or iteration number (``iterPeriod``).
-For instance, to print kinetic energy (using :yref:`kineticEnergy<yade._utils.kineticEnergy>`) every 5 seconds, this engine will be put to ``O.engines``::
+For instance, to print kinetic energy (using :yref:`kineticEnergy<yade.utils.kineticEnergy>`) every 5 seconds, the following engine will be put to ``O.engines``::
+
PyRunner(command="print 'kinetic energy',kineticEnergy()",realPeriod=5)
For running more complex commands, it is convenient to define an external function and only call it from within the engine. Since the ``command`` is run in the script's namespace, functions defined within scripts can be called. Let us print information on interaction between bodies 0 and 1 periodically::
@@ -1138,7 +1139,7 @@
Frequently, decisions have to be made based on evolution of the simulation itself, which is not yet known. In such case, a function checking some specific condition is called periodically; if the condition is satisfied, ``O.pause`` or other functions can be called to stop the stimulation. See documentation for :yref:`Omega.run`, :yref:`Omega.pause`, :yref:`Omega.step`, :yref:`Omega.stopAtIter` for details.
-For simulations that seek static equilibrium, the :yref:`unbalancedForce<yade._utils.unbalancedForce>` can provide a useful metrics (see its documentation for details); for a desired value of ``1e-2`` or less, for instance, we can use::
+For simulations that seek static equilibrium, the :yref:`unbalancedForce<yade.utils.unbalancedForce>` can provide a useful metrics (see its documentation for details); for a desired value of ``1e-2`` or less, for instance, we can use::
def checkUnbalanced():
@@ -1557,7 +1558,7 @@
Micro-stress
------------
-Stress fields can be generated by combining the volume returned by TesselationWrapper to per-particle stress given by :yref:`bodyStressTensors<yade._utils.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 [Catalano2014a]_ (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:
+Stress fields can be generated by combining the volume returned by TesselationWrapper to per-particle stress given by :yref:`bodyStressTensors<yade.utils.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 [Catalano2014a]_ (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 'gui/qt4/GLViewer.cpp'
--- gui/qt4/GLViewer.cpp 2014-07-03 17:20:40 +0000
+++ gui/qt4/GLViewer.cpp 2014-07-14 20:08:34 +0000
@@ -16,7 +16,6 @@
#include<yade/core/Scene.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/DisplayParameters.hpp>
-#include<boost/filesystem/operations.hpp>
#include<boost/algorithm/string.hpp>
#include<sstream>
#include<iomanip>
=== modified file 'gui/qt4/GLViewer.hpp'
--- gui/qt4/GLViewer.hpp 2014-07-03 07:15:45 +0000
+++ gui/qt4/GLViewer.hpp 2014-07-14 20:08:34 +0000
@@ -12,7 +12,6 @@
#include<QGLViewer/qglviewer.h>
#include<QGLViewer/manipulatedFrame.h>
#include<QGLViewer/constraint.h>
-#include<set>
using std::setw;
using std::setfill;
=== modified file 'gui/qt4/GLViewerDisplay.cpp'
--- gui/qt4/GLViewerDisplay.cpp 2014-07-03 17:24:54 +0000
+++ gui/qt4/GLViewerDisplay.cpp 2014-07-14 20:08:34 +0000
@@ -16,7 +16,6 @@
#include<yade/core/Scene.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/DisplayParameters.hpp>
-#include<boost/filesystem/operations.hpp>
#include<boost/algorithm/string.hpp>
#include<sstream>
#include<iomanip>
@@ -40,7 +39,7 @@
const shared_ptr<DisplayParameters>& dp=dispParams[n];
string val;
if(dp->getValue("OpenGLRenderer",val)){ istringstream oglre(val);
- yade::ObjectIO::load<TYPEOF(renderer),boost::archive::xml_iarchive>(oglre,"renderer",renderer);
+ yade::ObjectIO::load<decltype(renderer),boost::archive::xml_iarchive>(oglre,"renderer",renderer);
}
else { LOG_WARN("OpenGLRenderer configuration not found in display parameters, skipped.");}
if(dp->getValue("GLViewer",val)){ GLViewer::setState(val); displayMessage("Loaded view configuration #"+boost::lexical_cast<string>(n)); }
@@ -53,7 +52,7 @@
if(dispParams.size()<=n){while(dispParams.size()<=n) dispParams.push_back(shared_ptr<DisplayParameters>(new DisplayParameters));} assert(n<dispParams.size());
shared_ptr<DisplayParameters>& dp=dispParams[n];
ostringstream oglre;
- yade::ObjectIO::save<TYPEOF(renderer),boost::archive::xml_oarchive>(oglre,"renderer",renderer);
+ yade::ObjectIO::save<decltype(renderer),boost::archive::xml_oarchive>(oglre,"renderer",renderer);
dp->setValue("OpenGLRenderer",oglre.str());
dp->setValue("GLViewer",GLViewer::getState());
displayMessage("Saved view configuration ot #"+boost::lexical_cast<string>(n));
=== modified file 'gui/qt4/GLViewerMouse.cpp'
--- gui/qt4/GLViewerMouse.cpp 2014-07-03 17:20:40 +0000
+++ gui/qt4/GLViewerMouse.cpp 2014-07-14 20:08:34 +0000
@@ -16,7 +16,6 @@
#include<yade/core/Scene.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/DisplayParameters.hpp>
-#include<boost/filesystem/operations.hpp>
#include<boost/algorithm/string.hpp>
#include<sstream>
#include<iomanip>
=== modified file 'lib/base/Math.hpp'
--- lib/base/Math.hpp 2014-07-04 17:46:32 +0000
+++ lib/base/Math.hpp 2014-07-14 20:08:34 +0000
@@ -63,6 +63,11 @@
#include <boost/python/class.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/foreach.hpp>
+#include <boost/tuple/tuple.hpp>
+#include <boost/filesystem/convenience.hpp>
+#include <boost/filesystem/operations.hpp>
+#include <boost/filesystem/exception.hpp>
+#include <boost/numeric/conversion/bounds.hpp>
using boost::shared_ptr;
@@ -70,6 +75,18 @@
#define FOREACH BOOST_FOREACH
#endif
+#ifndef YADE_PTR_CAST
+ #define YADE_PTR_CAST boost::static_pointer_cast
+#endif
+
+#ifndef YADE_CAST
+ #define YADE_CAST static_cast
+#endif
+
+#ifndef YADE_PTR_DYN_CAST
+ #define YADE_PTR_DYN_CAST boost::dynamic_pointer_cast
+#endif
+
#define EIGEN_DONT_PARALLELIZE
#ifdef YADE_MASK_ARBITRARY
=== modified file 'lib/factory/DynLibManager.cpp'
--- lib/factory/DynLibManager.cpp 2014-07-03 17:20:40 +0000
+++ lib/factory/DynLibManager.cpp 2014-07-14 20:08:34 +0000
@@ -8,10 +8,6 @@
*************************************************************************/
#include "DynLibManager.hpp"
-
-#include<boost/filesystem/operations.hpp>
-#include<boost/filesystem/convenience.hpp>
-
#include "ClassFactory.hpp"
CREATE_LOGGER(DynLibManager);
=== modified file 'lib/serialization/Serializable.hpp'
--- lib/serialization/Serializable.hpp 2014-07-03 17:24:54 +0000
+++ lib/serialization/Serializable.hpp 2014-07-11 17:39:21 +0000
@@ -75,15 +75,15 @@
#define _DEF_READWRITE_BY_VALUE(thisClass,attr,doc) add_property(/*attr name*/BOOST_PP_STRINGIZE(attr),/*read access*/boost::python::make_getter(&thisClass::attr,boost::python::return_value_policy<boost::python::return_by_value>()),/*write access*/boost::python::make_setter(&thisClass::attr,boost::python::return_value_policy<boost::python::return_by_value>()),/*docstring*/doc)
// not sure if this is correct: the getter works by value, the setter by reference (the default)...?
-#define _DEF_READWRITE_BY_VALUE_POSTLOAD(thisClass,attr,doc) add_property(/*attr name*/BOOST_PP_STRINGIZE(attr),/*read access*/boost::python::make_getter(&thisClass::attr,boost::python::return_value_policy<boost::python::return_by_value>()),/*write access*/ make_setter_postLoad<thisClass,TYPEOF(thisClass::attr),&thisClass::attr>,/*docstring*/doc)
+#define _DEF_READWRITE_BY_VALUE_POSTLOAD(thisClass,attr,doc) add_property(/*attr name*/BOOST_PP_STRINGIZE(attr),/*read access*/boost::python::make_getter(&thisClass::attr,boost::python::return_value_policy<boost::python::return_by_value>()),/*write access*/ make_setter_postLoad<thisClass,decltype(thisClass::attr),&thisClass::attr>,/*docstring*/doc)
#define _DEF_READONLY_BY_VALUE(thisClass,attr,doc) add_property(/*attr name*/BOOST_PP_STRINGIZE(attr),/*read access*/boost::python::make_getter(&thisClass::attr,boost::python::return_value_policy<boost::python::return_by_value>()),/*docstring*/doc)
/* Huh, add_static_property does not support doc argument (add_property does); if so, use add_property for now at least... */
#define _DEF_READWRITE_BY_VALUE_STATIC(thisClass,attr,doc) _DEF_READWRITE_BY_VALUE(thisClass,attr,doc)
// the conditional yade::py_wrap_ref should be eliminated by compiler at compile-time, as it depends only on types, not their values
// most of this could be written with templates, including flags (ints can be template args)
-#define _DEF_READWRITE_CUSTOM(thisClass,attr) if(!(_ATTR_FLG(attr) & yade::Attr::hidden)){ bool _ro(_ATTR_FLG(attr) & Attr::readonly), _post(_ATTR_FLG(attr) & Attr::triggerPostLoad), _ref(yade::py_wrap_ref<TYPEOF(thisClass::_ATTR_NAM(attr))>::value); std::string docStr(_ATTR_DOC(attr)); docStr+=" :yattrflags:`"+boost::lexical_cast<string>(_ATTR_FLG(attr))+"` "; \
+#define _DEF_READWRITE_CUSTOM(thisClass,attr) if(!(_ATTR_FLG(attr) & yade::Attr::hidden)){ bool _ro(_ATTR_FLG(attr) & Attr::readonly), _post(_ATTR_FLG(attr) & Attr::triggerPostLoad), _ref(yade::py_wrap_ref<decltype(thisClass::_ATTR_NAM(attr))>::value); std::string docStr(_ATTR_DOC(attr)); docStr+=" :yattrflags:`"+boost::lexical_cast<string>(_ATTR_FLG(attr))+"` "; \
if ( _ref && !_ro && !_post) _classObj.def_readwrite(_ATTR_NAM_STR(attr),&thisClass::_ATTR_NAM(attr),docStr.c_str()); \
- else if ( _ref && !_ro && _post) _classObj.add_property(_ATTR_NAM_STR(attr),boost::python::make_getter(&thisClass::_ATTR_NAM(attr)),make_setter_postLoad<thisClass,TYPEOF(thisClass::_ATTR_NAM(attr)),&thisClass::_ATTR_NAM(attr)>,docStr.c_str()); \
+ else if ( _ref && !_ro && _post) _classObj.add_property(_ATTR_NAM_STR(attr),boost::python::make_getter(&thisClass::_ATTR_NAM(attr)),make_setter_postLoad<thisClass,decltype(thisClass::_ATTR_NAM(attr)),&thisClass::_ATTR_NAM(attr)>,docStr.c_str()); \
else if ( _ref && _ro) _classObj.def_readonly(_ATTR_NAM_STR(attr),&thisClass::_ATTR_NAM(attr),docStr.c_str()); \
else if (!_ref && !_ro && !_post) _classObj._DEF_READWRITE_BY_VALUE(thisClass,_ATTR_NAM(attr),docStr.c_str()); \
else if (!_ref && !_ro && _post) _classObj._DEF_READWRITE_BY_VALUE_POSTLOAD(thisClass,_ATTR_NAM(attr),docStr.c_str()); \
@@ -98,11 +98,11 @@
// gcc<=4.3 is not able to compile this code; we will just not generate any code for deprecated attributes in such case
#if !defined(__GNUG__) || (defined(__GNUG__) && (__GNUC__ > 4 || (__GNUC__==4 && __GNUC_MINOR__ > 3)))
// gcc > 4.3 && non-gcc compilers
- #define _PYSET_ATTR_DEPREC(x,thisClass,z) if(key==BOOST_PP_STRINGIZE(_DEPREC_OLDNAME(z))){ _DEPREC_WARN(thisClass,z); _DEPREC_NEWNAME(z)=boost::python::extract<TYPEOF(_DEPREC_NEWNAME(z))>(value); return; }
+ #define _PYSET_ATTR_DEPREC(x,thisClass,z) if(key==BOOST_PP_STRINGIZE(_DEPREC_OLDNAME(z))){ _DEPREC_WARN(thisClass,z); _DEPREC_NEWNAME(z)=boost::python::extract<decltype(_DEPREC_NEWNAME(z))>(value); return; }
#define _PYATTR_DEPREC_DEF(x,thisClass,z) .add_property(BOOST_PP_STRINGIZE(_DEPREC_OLDNAME(z)),&thisClass::BOOST_PP_CAT(_getDeprec_,_DEPREC_OLDNAME(z)),&thisClass::BOOST_PP_CAT(_setDeprec_,_DEPREC_OLDNAME(z)),"|ydeprecated| alias for :yref:`" BOOST_PP_STRINGIZE(_DEPREC_NEWNAME(z)) "<" BOOST_PP_STRINGIZE(thisClass) "." BOOST_PP_STRINGIZE(_DEPREC_NEWNAME(z)) ">` (" _DEPREC_COMMENT(z) ")")
#define _PYHASKEY_ATTR_DEPREC(x,thisClass,z) if(key==BOOST_PP_STRINGIZE(_DEPREC_OLDNAME(z))) return true;
/* accessors functions ussing warning */
- #define _ACCESS_DEPREC(x,thisClass,z) /*getter*/ TYPEOF(_DEPREC_NEWNAME(z)) BOOST_PP_CAT(_getDeprec_,_DEPREC_OLDNAME(z))(){_DEPREC_WARN(thisClass,z); return _DEPREC_NEWNAME(z); } /*setter*/ void BOOST_PP_CAT(_setDeprec_,_DEPREC_OLDNAME(z))(const TYPEOF(_DEPREC_NEWNAME(z))& val){_DEPREC_WARN(thisClass,z); _DEPREC_NEWNAME(z)=val; }
+ #define _ACCESS_DEPREC(x,thisClass,z) /*getter*/ decltype(_DEPREC_NEWNAME(z)) BOOST_PP_CAT(_getDeprec_,_DEPREC_OLDNAME(z))(){_DEPREC_WARN(thisClass,z); return _DEPREC_NEWNAME(z); } /*setter*/ void BOOST_PP_CAT(_setDeprec_,_DEPREC_OLDNAME(z))(const decltype(_DEPREC_NEWNAME(z))& val){_DEPREC_WARN(thisClass,z); _DEPREC_NEWNAME(z)=val; }
#else
#define _PYSET_ATTR_DEPREC(x,y,z)
#define _PYATTR_DEPREC_DEF(x,y,z)
@@ -114,7 +114,7 @@
// loop bodies for attribute access
#define _PYGET_ATTR(x,y,z) if(key==_ATTR_NAM_STR(z)) return boost::python::object(_ATTR_NAM(z));
//#define _PYSET_ATTR(x,y,z) if(key==_ATTR_NAM_STR(z)) { _ATTR_NAM(z)=boost::python::extract<typeof(_ATTR_NAM(z))>(t[1]); boost::python::delitem(d,boost::python::object(_ATTR_NAM(z))); continue; }
-#define _PYSET_ATTR(x,y,z) if(key==_ATTR_NAM_STR(z)) { _ATTR_NAM(z)=boost::python::extract<TYPEOF(_ATTR_NAM(z))>(value); return; }
+#define _PYSET_ATTR(x,y,z) if(key==_ATTR_NAM_STR(z)) { _ATTR_NAM(z)=boost::python::extract<decltype(_ATTR_NAM(z))>(value); return; }
#define _PYKEYS_ATTR(x,y,z) ret.append(_ATTR_NAM_STR(z));
#define _PYHASKEY_ATTR(x,y,z) if(key==_ATTR_NAM_STR(z)) return true;
#define _PYDICT_ATTR(x,y,z) if(!(_ATTR_FLG(z) & yade::Attr::hidden)) ret[_ATTR_NAM_STR(z)]=boost::python::object(_ATTR_NAM(z));
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-07-02 16:18:24 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-07-16 16:49:41 +0000
@@ -116,7 +116,7 @@
//This is the influx term
if (cell->info().Pcondition) cell->info().averageVelocity() = cell->info().averageVelocity() - (totFlowRate)*((Point) cell->info()-CGAL::ORIGIN );
//now divide by volume
- cell->info().averageVelocity() = cell->info().averageVelocity() /abs(cell->info().volume());
+ cell->info().averageVelocity() = cell->info().averageVelocity() /std::abs(cell->info().volume());
}
}
@@ -239,7 +239,7 @@
double pressure = 0.f;
int cell=0;
for (int i=0; i<captures; i++){
- for (double Z=min(zMin,zMax); Z<=max(zMin,zMax); Z+=abs(Rz)) {
+ for (double Z=min(zMin,zMax); Z<=max(zMin,zMax); Z+=std::abs(Rz)) {
permeameter = Tri.locate(Point(X, Y, Z));
pressure+=permeameter->info().p();
cell++;
@@ -318,7 +318,7 @@
CVector fluidSurfk = cell->info().facetSurfaces[j]*cell->info().facetFluidSurfacesRatio[j];
/// handle fictious vertex since we can get the projected surface easily here
if (cell->vertex(j)->info().isFictious) {
- Real projSurf=abs(Surfk[boundary(cell->vertex(j)->info().id()).coordinate]);
+ Real projSurf=std::abs(Surfk[boundary(cell->vertex(j)->info().id()).coordinate]);
tempVect=-projSurf*boundary(cell->vertex(j)->info().id()).normal;
cell->vertex(j)->info().forces = cell->vertex(j)->info().forces+tempVect*cell->info().p();
//define the cached value for later use with cache*p
@@ -526,7 +526,7 @@
if (S0==0) S0=checkSphereFacetOverlap(v1,v2,v0);
if (S0==0) S0=checkSphereFacetOverlap(v2,v0,v1);
//take absolute value, since in rare cases the surface can be negative (overlaping spheres)
- fluidArea=abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
+ fluidArea=std::abs(area-crossSections[0]-crossSections[1]-crossSections[2]+S0);
cell->info().facetFluidSurfacesRatio[j]=fluidArea/area;
k=(fluidArea * pow(radius,2)) / (8*viscosity*distance);
meanDistance += distance;
@@ -988,7 +988,7 @@
double viscosity = viscosity;
double gravity = 1;
double Vdarcy = Q1/Section;
- double DeltaP = abs(PInf-PSup);
+ double DeltaP = std::abs(PInf-PSup);
double DeltaH = DeltaP/ (density*gravity);
double k = viscosity*Vdarcy*DeltaY / DeltaP; /**m²**/
double Ks = k*(density*gravity)/viscosity; /**m/s**/
@@ -1185,7 +1185,7 @@
boundary(yMaxId).flowCondition=0;
boundary(yMinId).value=0;
boundary(yMaxId).value=1;
- double pZero = abs((boundary(yMinId).value-boundary(yMaxId).value)/2);
+ double pZero = std::abs((boundary(yMinId).value-boundary(yMaxId).value)/2);
initializePressure( pZero );
gaussSeidel();
const char *kk = "Permeability";
@@ -1214,8 +1214,8 @@
double Ry = (yMax-yMin) /intervals;
double Rz = (zMax-zMin) /intervals;
double X=0.5;
- for (double Y=min(yMax,yMin); Y<=max(yMax,yMin); Y=Y+abs(Ry)) {
- for (double Z=min(zMin,zMax); Z<=max(zMin,zMax); Z=Z+abs(Rz)) {
+ for (double Y=min(yMax,yMin); Y<=max(yMax,yMin); Y=Y+std::abs(Ry)) {
+ for (double Z=min(zMin,zMax); Z<=max(zMin,zMax); Z=Z+std::abs(Rz)) {
permeameter = Tri.locate(Point(X, Y, Z));
consFile << permeameter->info().p() <<" ";
}
=== modified file 'lib/triangulation/Network.ipp'
--- lib/triangulation/Network.ipp 2014-07-09 07:07:14 +0000
+++ lib/triangulation/Network.ipp 2014-07-17 04:21:07 +0000
@@ -15,6 +15,8 @@
#define FAST
namespace CGT {
+
+using std::abs;
// template<class Tesselation> const double Network<Tesselation>::FAR = 50000;
template<class Tesselation> const double Network<Tesselation>::ONE_THIRD = 1.0/3.0;
=== modified file 'pkg/common/Cylinder.hpp'
--- pkg/common/Cylinder.hpp 2014-05-23 13:05:19 +0000
+++ pkg/common/Cylinder.hpp 2014-07-14 20:08:33 +0000
@@ -10,9 +10,7 @@
#ifdef YADE_OPENGL
#include<yade/pkg/common/GLDrawFunctors.hpp>
#endif
-#include<yade/pkg/dem/CohFrictPhys.hpp>
-#include<yade/pkg/dem/CohFrictMat.hpp>
-#include<yade/pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp>
+#include<yade/pkg/dem/CohesiveFrictionalContactLaw.hpp>
#include<yade/pkg/common/CylScGeom6D.hpp>
class Cylinder: public Sphere{
=== modified file 'pkg/common/Grid.hpp'
--- pkg/common/Grid.hpp 2014-05-23 13:05:19 +0000
+++ pkg/common/Grid.hpp 2014-07-14 20:08:33 +0000
@@ -23,7 +23,7 @@
#pragma once
#include "Sphere.hpp"
#include <yade/pkg/dem/FrictPhys.hpp>
-#include <yade/pkg/dem/CohFrictPhys.hpp>
+#include <yade/pkg/dem/CohesiveFrictionalContactLaw.hpp>
#include <yade/pkg/dem/ScGeom.hpp>
#include <yade/core/Body.hpp>
#include <yade/pkg/common/Dispatching.hpp>
=== modified file 'pkg/common/PersistentTriangulationCollider.hpp'
--- pkg/common/PersistentTriangulationCollider.hpp 2010-11-07 11:46:20 +0000
+++ pkg/common/PersistentTriangulationCollider.hpp 2014-07-14 20:08:34 +0000
@@ -9,10 +9,6 @@
#pragma once
#include<yade/pkg/common/Collider.hpp>
#include<yade/core/InteractionContainer.hpp>
-#include <list>
-#include <set>
-#include <vector>
-#include <algorithm>
#include<yade/pkg/dem/TesselationWrapper.hpp>
/*! \brief Collision detection engine based on regular triangulation.
=== modified file 'pkg/dem/CapillaryPhys.cpp'
--- pkg/dem/CapillaryPhys.cpp 2010-11-12 08:03:16 +0000
+++ pkg/dem/CapillaryPhys.cpp 2014-07-14 20:08:34 +0000
@@ -1,8 +1,39 @@
#include <yade/pkg/dem/CapillaryPhys.hpp>
-
-CapillaryPhys::~CapillaryPhys()
+#include <yade/pkg/dem/ScGeom.hpp>
+
+YADE_PLUGIN((CapillaryPhys)(Ip2_FrictMat_FrictMat_CapillaryPhys));
+
+// The following code was moved from Ip2_FrictMat_FrictMat_CapillaryPhys.cpp
+
+void Ip2_FrictMat_FrictMat_CapillaryPhys::go( const shared_ptr<Material>& b1 //FrictMat
+ , const shared_ptr<Material>& b2 // FrictMat
+ , const shared_ptr<Interaction>& interaction)
{
-}
-
-YADE_PLUGIN((CapillaryPhys));
-
+ ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
+ if(geom)
+ {
+ if(!interaction->phys)
+ {
+ const shared_ptr<FrictMat>& sdec1 = YADE_PTR_CAST<FrictMat>(b1);
+ const shared_ptr<FrictMat>& sdec2 = YADE_PTR_CAST<FrictMat>(b2);
+
+ if (!interaction->phys) interaction->phys = shared_ptr<CapillaryPhys>(new CapillaryPhys());
+ const shared_ptr<CapillaryPhys>& contactPhysics = YADE_PTR_CAST<CapillaryPhys>(interaction->phys);
+
+ Real Ea = sdec1->young;
+ Real Eb = sdec2->young;
+ Real Va = sdec1->poisson;
+ Real Vb = sdec2->poisson;
+ Real Da = geom->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size)
+ Real Db = geom->radius2; // FIXME - as above
+ Real fa = sdec1->frictionAngle;
+ Real fb = sdec2->frictionAngle;
+ Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
+ Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
+
+ contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
+ }
+ }
+};
=== modified file 'pkg/dem/CapillaryPhys.hpp'
--- pkg/dem/CapillaryPhys.hpp 2013-05-30 20:17:45 +0000
+++ pkg/dem/CapillaryPhys.hpp 2014-07-14 20:08:34 +0000
@@ -7,13 +7,15 @@
*************************************************************************/
#pragma once
#include<yade/pkg/dem/FrictPhys.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
class CapillaryPhys : public FrictPhys
{
public :
int currentIndexes [4]; // used for faster interpolation (stores previous positions in tables)
- virtual ~CapillaryPhys();
+ virtual ~CapillaryPhys() {};
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(CapillaryPhys,FrictPhys,"Physics (of interaction) for Law2_ScGeom_CapillaryPhys_Capillarity.",
((bool,meniscus,false,,"Presence of a meniscus if true"))
@@ -36,5 +38,16 @@
REGISTER_SERIALIZABLE(CapillaryPhys);
+class Ip2_FrictMat_FrictMat_CapillaryPhys : public IPhysFunctor
+{
+ public :
+ virtual void go( const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ FUNCTOR2D(FrictMat,FrictMat);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_FrictMat_CapillaryPhys,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",,;);
+
+};
+REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_CapillaryPhys);
=== modified file 'pkg/dem/CapillaryTriaxialTest.cpp'
--- pkg/dem/CapillaryTriaxialTest.cpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/CapillaryTriaxialTest.cpp 2014-07-14 20:08:34 +0000
@@ -11,8 +11,7 @@
#include<yade/pkg/dem/ElasticContactLaw.hpp>
#include <yade/pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp>
-// #include<yade/pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp>
-#include<yade/pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.hpp>
+#include<yade/pkg/dem/CapillaryPhys.hpp>
#include<yade/pkg/common/ElastMat.hpp>
#include<yade/pkg/dem/GlobalStiffnessTimeStepper.hpp>
@@ -42,7 +41,6 @@
#include<yade/pkg/dem/Shop.hpp>
-#include <boost/filesystem/convenience.hpp>
#include <boost/numeric/conversion/bounds.hpp>
#include <boost/limits.hpp>
=== removed file 'pkg/dem/CohFrictMat.cpp'
--- pkg/dem/CohFrictMat.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/CohFrictMat.cpp 1970-01-01 00:00:00 +0000
@@ -1,19 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include "CohFrictMat.hpp"
-
-
-CohFrictMat::~CohFrictMat()
-{
-}
-
-YADE_PLUGIN((CohFrictMat));
-
-
-
=== removed file 'pkg/dem/CohFrictMat.hpp'
--- pkg/dem/CohFrictMat.hpp 2014-06-02 09:07:27 +0000
+++ pkg/dem/CohFrictMat.hpp 1970-01-01 00:00:00 +0000
@@ -1,40 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#pragma once
-
-
-#include<yade/pkg/common/ElastMat.hpp>
-
-
-class CohFrictMat : public FrictMat
-{
- public :
- virtual ~CohFrictMat ();
-
-/// Serialization
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictMat,FrictMat,"",
- ((bool,isCohesive,true,,""))
- ((Real,alphaKr,2.0,,"Dimensionless rolling stiffness."))
- ((Real,alphaKtw,2.0,,"Dimensionless twist stiffness."))
- ((Real,etaRoll,-1.,,"Dimensionless rolling (aka 'bending') strength. If negative, rolling moment will be elastic."))
- ((Real,etaTwist,-1.,,"Dimensionless twisting strength. If negative, twist moment will be elastic."))
- ((Real,normalCohesion,-1,,"Tensile strength, homogeneous to a pressure. If negative the normal force is purely elastic."))
- ((Real,shearCohesion,-1,,"Shear strength, homogeneous to a pressure. If negative the shear force is purely elastic."))
- ((bool,momentRotationLaw,false,,"Use bending/twisting moment at contact. The contact will have moments only if both bodies have this flag true. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))
- ,
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(CohFrictMat,FrictMat);
-};
-
-REGISTER_SERIALIZABLE(CohFrictMat);
-
-
=== removed file 'pkg/dem/CohFrictPhys.cpp'
--- pkg/dem/CohFrictPhys.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/CohFrictPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,23 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include "CohFrictPhys.hpp"
-
-void CohFrictPhys::SetBreakingState()
-{
- cohesionBroken = true;
- normalAdhesion = 0;
- shearAdhesion = 0;
-}
-
-CohFrictPhys::~CohFrictPhys()
-{
-}
-YADE_PLUGIN((CohFrictPhys));
-
-
=== removed file 'pkg/dem/CohFrictPhys.hpp'
--- pkg/dem/CohFrictPhys.hpp 2014-05-26 16:49:23 +0000
+++ pkg/dem/CohFrictPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,47 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/common/NormShearPhys.hpp>
-#include<yade/pkg/dem/FrictPhys.hpp>
-
-class CohFrictPhys : public FrictPhys
-{
- public :
- virtual ~CohFrictPhys();
- void SetBreakingState();
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictPhys,FrictPhys,"",
- ((bool,cohesionDisablesFriction,false,,"is shear strength the sum of friction and adhesion or only adhesion?"))
- ((bool,cohesionBroken,true,,"is cohesion active? Set to false at the creation of a cohesive contact, and set to true when a fragile contact is broken"))
- ((bool,fragile,true,,"do cohesion disapear when contact strength is exceeded?"))
- ((Real,kr,0,,"rotational stiffness [N.m/rad]"))
- ((Real,ktw,0,,"twist stiffness [N.m/rad]"))
- ((Real,maxRollPl,0.0,,"Coefficient of rolling friction (negative means elastic)."))
- ((Real,maxTwistPl,0.0,,"Coefficient of twisting friction (negative means elastic)."))
- ((Real,normalAdhesion,0,,"tensile strength"))
- ((Real,shearAdhesion,0,,"cohesive part of the shear strength (a frictional term might be added depending on :yref:`CohFrictPhys::cohesionDisablesFriction`)"))
- ((Real,unp,0,,"plastic normal displacement, only used for tensile behaviour and if :yref:`CohFrictPhys::fragile` =false."))
- ((Real,unpMax,0,,"maximum value of plastic normal displacement, after that the interaction breaks even if :yref:`CohFrictPhys::fragile` =false. The default value (0) means no maximum."))
- ((bool,momentRotationLaw,false,,"use bending/twisting moment at contacts. See :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::always_use_moment_law` for details."))
- ((bool,initCohesion,false,,"Initialize the cohesive behaviour with current state as equilibrium state (same as :yref:`Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::setCohesionNow` but acting on only one interaction)"))
- ((Real,creep_viscosity,-1,,"creep viscosity [Pa.s/m]."))
- // internal attributes
- ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment"))
- ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment"))
- ,
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(CohFrictPhys,FrictPhys);
-
-};
-
-REGISTER_SERIALIZABLE(CohFrictPhys);
-
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp 2014-06-02 09:07:27 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp 2014-07-14 20:08:33 +0000
@@ -7,13 +7,11 @@
*************************************************************************/
#include "CohesiveFrictionalContactLaw.hpp"
-#include<yade/pkg/dem/CohFrictMat.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/CohFrictPhys.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-YADE_PLUGIN((CohesiveFrictionalContactLaw)(Law2_ScGeom6D_CohFrictPhys_CohesionMoment));
+YADE_PLUGIN((CohesiveFrictionalContactLaw)(Law2_ScGeom6D_CohFrictPhys_CohesionMoment)(CohFrictMat)(CohFrictPhys)(Ip2_CohFrictMat_CohFrictMat_CohFrictPhys));
CREATE_LOGGER(Law2_ScGeom6D_CohFrictPhys_CohesionMoment);
Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::getPlasticDissipation() {return (Real) plasticDissipation;}
@@ -226,3 +224,76 @@
/// Moment law END ///
}
}
+
+
+void Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::go(const shared_ptr<Material>& b1 // CohFrictMat
+ , const shared_ptr<Material>& b2 // CohFrictMat
+ , const shared_ptr<Interaction>& interaction)
+{
+ CohFrictMat* sdec1 = static_cast<CohFrictMat*>(b1.get());
+ CohFrictMat* sdec2 = static_cast<CohFrictMat*>(b2.get());
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
+
+ //Create cohesive interractions only once
+ if (setCohesionNow && cohesionDefinitionIteration==-1) cohesionDefinitionIteration=scene->iter;
+ if (setCohesionNow && cohesionDefinitionIteration!=-1 && cohesionDefinitionIteration!=scene->iter) {
+ cohesionDefinitionIteration = -1;
+ setCohesionNow = 0;}
+
+ if (geom) {
+ if (!interaction->phys) {
+ interaction->phys = shared_ptr<CohFrictPhys>(new CohFrictPhys());
+ CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());
+ Real Ea = sdec1->young;
+ Real Eb = sdec2->young;
+ Real Va = sdec1->poisson;
+ Real Vb = sdec2->poisson;
+ Real Da = geom->radius1;
+ Real Db = geom->radius2;
+ Real fa = sdec1->frictionAngle;
+ Real fb = sdec2->frictionAngle;
+ Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
+
+ // harmonic average of alphas parameters
+ Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
+ Real AlphaKtw;
+ if (sdec1->alphaKtw && sdec2->alphaKtw) AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
+ else AlphaKtw=0;
+
+ Real Ks;
+ if (Va && Vb) Ks = 2.0*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
+ else Ks=0;
+
+ contactPhysics->kr = Da*Db*Ks*AlphaKr;
+ contactPhysics->ktw = Da*Db*Ks*AlphaKtw;
+ contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
+
+ if ((setCohesionOnNewContacts || setCohesionNow) && sdec1->isCohesive && sdec2->isCohesive)
+ {
+ contactPhysics->cohesionBroken = false;
+ contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(Db, Da),2);
+ contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(Db, Da),2);
+ geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
+ }
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
+
+ contactPhysics->maxRollPl = min(sdec1->etaRoll*Da,sdec2->etaRoll*Db);
+ contactPhysics->maxTwistPl = min(sdec1->etaTwist*Da,sdec2->etaTwist*Db);
+ contactPhysics->momentRotationLaw=(sdec1->momentRotationLaw && sdec2->momentRotationLaw);
+ }
+ else {// !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
+ CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());
+ if ((setCohesionNow && sdec1->isCohesive && sdec2->isCohesive) || contactPhysics->initCohesion)
+ {
+ contactPhysics->cohesionBroken = false;
+ contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
+ contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
+
+ geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
+ contactPhysics->initCohesion=false;
+ }
+ }
+ }
+};
+
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-07-14 20:08:33 +0000
@@ -9,10 +9,72 @@
#pragma once
#include<yade/core/GlobalEngine.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/CohFrictPhys.hpp>
-#include<yade/pkg/common/Dispatching.hpp>
#include<boost/tuple/tuple.hpp>
+#include<yade/pkg/common/NormShearPhys.hpp>
+#include<yade/pkg/dem/FrictPhys.hpp>
+
+
+// The following code was moved from CohFrictMat.hpp
+class CohFrictMat : public FrictMat
+{
+ public :
+ virtual ~CohFrictMat () {};
+
+/// Serialization
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictMat,FrictMat,"",
+ ((bool,isCohesive,true,,""))
+ ((Real,alphaKr,2.0,,"Dimensionless rolling stiffness."))
+ ((Real,alphaKtw,2.0,,"Dimensionless twist stiffness."))
+ ((Real,etaRoll,-1.,,"Dimensionless rolling (aka 'bending') strength. If negative, rolling moment will be elastic."))
+ ((Real,etaTwist,-1.,,"Dimensionless twisting strength. If negative, twist moment will be elastic."))
+ ((Real,normalCohesion,-1,,"Tensile strength, homogeneous to a pressure. If negative the normal force is purely elastic."))
+ ((Real,shearCohesion,-1,,"Shear strength, homogeneous to a pressure. If negative the shear force is purely elastic."))
+ ((bool,momentRotationLaw,false,,"Use bending/twisting moment at contact. The contact will have moments only if both bodies have this flag true. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))
+ ,
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(CohFrictMat,FrictMat);
+};
+
+REGISTER_SERIALIZABLE(CohFrictMat);
+
+// The following code was moved from CohFrictPhys.hpp
+class CohFrictPhys : public FrictPhys
+{
+ public :
+ virtual ~CohFrictPhys() {};
+ void SetBreakingState() {cohesionBroken = true; normalAdhesion = 0; shearAdhesion = 0;};
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictPhys,FrictPhys,"",
+ ((bool,cohesionDisablesFriction,false,,"is shear strength the sum of friction and adhesion or only adhesion?"))
+ ((bool,cohesionBroken,true,,"is cohesion active? Set to false at the creation of a cohesive contact, and set to true when a fragile contact is broken"))
+ ((bool,fragile,true,,"do cohesion disapear when contact strength is exceeded?"))
+ ((Real,kr,0,,"rotational stiffness [N.m/rad]"))
+ ((Real,ktw,0,,"twist stiffness [N.m/rad]"))
+ ((Real,maxRollPl,0.0,,"Coefficient of rolling friction (negative means elastic)."))
+ ((Real,maxTwistPl,0.0,,"Coefficient of twisting friction (negative means elastic)."))
+ ((Real,normalAdhesion,0,,"tensile strength"))
+ ((Real,shearAdhesion,0,,"cohesive part of the shear strength (a frictional term might be added depending on :yref:`CohFrictPhys::cohesionDisablesFriction`)"))
+ ((Real,unp,0,,"plastic normal displacement, only used for tensile behaviour and if :yref:`CohFrictPhys::fragile` =false."))
+ ((Real,unpMax,0,,"maximum value of plastic normal displacement, after that the interaction breaks even if :yref:`CohFrictPhys::fragile` =false. The default value (0) means no maximum."))
+ ((bool,momentRotationLaw,false,,"use bending/twisting moment at contacts. See :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::always_use_moment_law` for details."))
+ ((bool,initCohesion,false,,"Initialize the cohesive behaviour with current state as equilibrium state (same as :yref:`Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::setCohesionNow` but acting on only one interaction)"))
+ ((Real,creep_viscosity,-1,,"creep viscosity [Pa.s/m]."))
+ // internal attributes
+ ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment"))
+ ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment"))
+ ,
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(CohFrictPhys,FrictPhys);
+};
+
+REGISTER_SERIALIZABLE(CohFrictPhys);
class Law2_ScGeom6D_CohFrictPhys_CohesionMoment: public LawFunctor{
public:
@@ -68,3 +130,23 @@
REGISTER_SERIALIZABLE(CohesiveFrictionalContactLaw);
+// The following code was moved from Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp
+class Ip2_CohFrictMat_CohFrictMat_CohFrictPhys : public IPhysFunctor
+{
+ public :
+ virtual void go( const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ int cohesionDefinitionIteration;
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_CohFrictMat_CohFrictMat_CohFrictPhys,IPhysFunctor,
+ "Generates cohesive-frictional interactions with moments, used in the contact law :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`. The normal/shear stiffness and friction definitions are the same as in :yref:`Ip2_FrictMat_FrictMat_FrictPhys`, check the documentation there for details.\n\nAdhesions related to the normal and the shear components are calculated form :yref:`CohFrictMat::normalCohesion` ($C_n$) and :yref:`CohFrictMat::shearlCohesion` ($C_s$). For particles of size $R_1$,$R_2$ the adhesion will be $a_i=C_i min(R_1,R_2)^2$, $i=n\\,s$.\n\nTwist and rolling stiffnesses are proportional to the shear stiffness through dimensionless factors alphaKtw and alphaKr, such that the rotational stiffnesses are defined by $k_s \\alpha_i R_1 R_2$, $i=tw\\,r$",
+ ((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts in current time-step. The flag is turned false automatically, so that assignment is done in the current timestep only."))
+ ((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts. If false, only existing contacts can be cohesive (also see :yref:`Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::setCohesionNow`), and new contacts are only frictional."))
+ ,
+ cohesionDefinitionIteration = -1;
+ );
+ FUNCTOR2D(CohFrictMat,CohFrictMat);
+};
+
+REGISTER_SERIALIZABLE(Ip2_CohFrictMat_CohFrictMat_CohFrictPhys);
=== modified file 'pkg/dem/CohesiveTriaxialTest.cpp'
--- pkg/dem/CohesiveTriaxialTest.cpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/CohesiveTriaxialTest.cpp 2014-07-14 20:08:34 +0000
@@ -9,8 +9,6 @@
#include "CohesiveTriaxialTest.hpp"
#include<yade/pkg/dem/CohesiveFrictionalContactLaw.hpp>
-#include<yade/pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp>
-#include<yade/pkg/dem/CohFrictMat.hpp>
#include<yade/pkg/dem/GlobalStiffnessTimeStepper.hpp>
#include<yade/pkg/dem/TriaxialStressController.hpp>
@@ -38,8 +36,6 @@
#include<yade/pkg/dem/Shop.hpp>
-#include <boost/filesystem/convenience.hpp>
-#include <boost/numeric/conversion/bounds.hpp>
#include <boost/limits.hpp>
// random
=== modified file 'pkg/dem/ConcretePM.hpp'
--- pkg/dem/ConcretePM.hpp 2014-04-10 10:25:02 +0000
+++ pkg/dem/ConcretePM.hpp 2014-07-14 20:08:34 +0000
@@ -53,7 +53,7 @@
#include<yade/pkg/common/NormShearPhys.hpp>
#include<yade/pkg/dem/DemXDofGeom.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp>
+#include<yade/pkg/dem/FrictPhys.hpp>
namespace py=boost::python;
=== modified file 'pkg/dem/ElasticContactLaw.hpp'
--- pkg/dem/ElasticContactLaw.hpp 2013-01-14 21:20:06 +0000
+++ pkg/dem/ElasticContactLaw.hpp 2014-07-14 20:08:34 +0000
@@ -13,9 +13,6 @@
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/lib/base/openmp-accu.hpp>
-#include<set>
-#include<boost/tuple/tuple.hpp>
-
class Law2_ScGeom_FrictPhys_CundallStrack: public LawFunctor{
public:
OpenMPAccumulator<Real> plasticDissipation;
=== modified file 'pkg/dem/FrictPhys.cpp'
--- pkg/dem/FrictPhys.cpp 2012-06-04 21:03:29 +0000
+++ pkg/dem/FrictPhys.cpp 2014-07-14 20:08:34 +0000
@@ -1,4 +1,69 @@
#include "FrictPhys.hpp"
-FrictPhys::~FrictPhys(){}
-ViscoFrictPhys::~ViscoFrictPhys(){}
-YADE_PLUGIN((FrictPhys)(ViscoFrictPhys));
+#include <yade/pkg/dem/ScGeom.hpp>
+YADE_PLUGIN((FrictPhys)(ViscoFrictPhys)(Ip2_FrictMat_FrictMat_ViscoFrictPhys)(Ip2_FrictMat_FrictMat_FrictPhys));
+
+// The following code was moved from Ip2_FrictMat_FrictMat_FrictPhys.hpp
+
+void Ip2_FrictMat_FrictMat_FrictPhys::go( const shared_ptr<Material>& b1
+ , const shared_ptr<Material>& b2
+ , const shared_ptr<Interaction>& interaction)
+{
+ if(interaction->phys) return;
+
+ const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
+ const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
+
+ Real Ra,Rb;//Vector3r normal;
+ assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode
+ GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
+ Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
+ Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
+
+ interaction->phys = shared_ptr<FrictPhys>(new FrictPhys());
+ const shared_ptr<FrictPhys>& contactPhysics = YADE_PTR_CAST<FrictPhys>(interaction->phys);
+ Real Ea = mat1->young;
+ Real Eb = mat2->young;
+ Real Va = mat1->poisson;
+ Real Vb = mat2->poisson;
+
+ //harmonic average of the two stiffnesses when (2*Ri*Ei) is the stiffness of a contact point on sphere "i"
+ Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
+ //same for shear stiffness
+ Real Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);
+
+ Real frictionAngle = (!frictAngle) ? std::min(mat1->frictionAngle,mat2->frictionAngle) : (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle);
+ contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
+};
+
+void Ip2_FrictMat_FrictMat_ViscoFrictPhys::go( const shared_ptr<Material>& b1
+ , const shared_ptr<Material>& b2
+ , const shared_ptr<Interaction>& interaction)
+{
+ if(interaction->phys) return;
+ const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
+ const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
+ interaction->phys = shared_ptr<ViscoFrictPhys>(new ViscoFrictPhys());
+ const shared_ptr<ViscoFrictPhys>& contactPhysics = YADE_PTR_CAST<ViscoFrictPhys>(interaction->phys);
+ Real Ea = mat1->young;
+ Real Eb = mat2->young;
+ Real Va = mat1->poisson;
+ Real Vb = mat2->poisson;
+
+ Real Ra,Rb;//Vector3r normal;
+ assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode
+ GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
+ Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
+ Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
+
+ //harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
+ Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
+ //same for shear stiffness
+ Real Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);
+
+ Real frictionAngle = (!frictAngle) ? std::min(mat1->frictionAngle,mat2->frictionAngle) : (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle);
+ contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
+};
=== modified file 'pkg/dem/FrictPhys.hpp'
--- pkg/dem/FrictPhys.hpp 2012-06-04 21:03:29 +0000
+++ pkg/dem/FrictPhys.hpp 2014-07-14 20:08:34 +0000
@@ -8,11 +8,14 @@
#pragma once
#include<yade/pkg/common/NormShearPhys.hpp>
+#include<yade/pkg/common/MatchMaker.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+#include<yade/pkg/common/Dispatching.hpp>
class FrictPhys: public NormShearPhys
{
public :
- virtual ~FrictPhys();
+ virtual ~FrictPhys() {};
YADE_CLASS_BASE_DOC_ATTRS_CTOR(FrictPhys,NormShearPhys,"The simple linear elastic-plastic interaction with friction angle, like in the traditional [CundallStrack1979]_",
((Real,tangensOfFrictionAngle,NaN,,"tan of angle of friction")),
createIndex()
@@ -24,7 +27,7 @@
class ViscoFrictPhys: public FrictPhys
{
public :
- virtual ~ViscoFrictPhys();
+ virtual ~ViscoFrictPhys() {};
YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscoFrictPhys,FrictPhys,"Temporary version of :yref:`FrictPhys` for compatibility with e.g. :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`",
((Vector3r,creepedShear,Vector3r(0,0,0),(Attr::readonly),"Creeped force (parallel)")),
createIndex()
@@ -32,3 +35,28 @@
REGISTER_CLASS_INDEX(ViscoFrictPhys,FrictPhys);
};
REGISTER_SERIALIZABLE(ViscoFrictPhys);
+
+// The following code was moved from Ip2_FrictMat_FrictMat_FrictPhys.hpp
+
+class Ip2_FrictMat_FrictMat_FrictPhys: public IPhysFunctor{
+ public:
+ virtual void go(const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ FUNCTOR2D(FrictMat,FrictMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_FrictPhys,IPhysFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under point load is defined here as $1/(E.D)$, with $E$ the stiffness of the sphere and $D$ its diameter. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. $1/(E_1.D_1)+1/(E_2.D_2)$ in the general case, or $2/(E.D)$ in the special case of equal sizes and equal stiffness. Note that summing compliances corresponds to an harmonic average of stiffnesss (as in e.g. [Scholtes2009a]_), which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor:\n\n $k_n = \\frac{E_1D_1*E_2D_2}{E_1D_1+E_2D_2}=\\frac{k_1*k_2}{k_1+k_2}$, with $k_i=E_iD_i$.\n\n The shear stiffness ks of one sphere is defined via the material parameter :yref:`ElastMat::poisson`, as ks=poisson*kn, and the resulting shear stiffness of the interaction will be also an harmonic average. In the case of a contact between a :yref:`ViscElMat` and a :yref:`FrictMat`, be sure to set :yref:`FrictMat::young` and :yref:`FrictMat::poisson`, otherwise the default value will be used.",
+ ((shared_ptr<MatchMaker>,frictAngle,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's friction angle. If ``None``, minimum value is used."))
+ );
+};
+REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_FrictPhys);
+
+class Ip2_FrictMat_FrictMat_ViscoFrictPhys: public Ip2_FrictMat_FrictMat_FrictPhys{
+ public:
+ virtual void go(const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ FUNCTOR2D(FrictMat,FrictMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_ViscoFrictPhys,Ip2_FrictMat_FrictMat_FrictPhys,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under symetric point loads is defined here as 1/(E.r), with E the stiffness of the sphere and r its radius, and corresponds to a compliance 1/(2.E.r)=1/(E.D) from each contact point. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. 1/(E.D1)+1/(E.D2) in the general case, or 1/(E.r) in the special case of equal sizes. Note that summing compliances corresponds to an harmonic average of stiffnesss, which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor.\n\nThe shear stiffness ks of one sphere is defined via the material parameter :yref:`ElastMat::poisson`, as ks=poisson*kn, and the resulting shear stiffness of the interaction will be also an harmonic average.",
+ );
+};
+REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_ViscoFrictPhys);
=== modified file 'pkg/dem/FrictViscoPM.hpp'
--- pkg/dem/FrictViscoPM.hpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/FrictViscoPM.hpp 2014-07-14 20:08:34 +0000
@@ -26,8 +26,6 @@
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/lib/base/openmp-accu.hpp>
-#include<boost/tuple/tuple.hpp>
-
/** This class holds information associated with each body */
class FrictViscoMat: public FrictMat {
public:
=== modified file 'pkg/dem/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp 2013-09-23 17:17:34 +0000
+++ pkg/dem/HertzMindlin.cpp 2014-07-14 20:08:33 +0000
@@ -1,18 +1,18 @@
// 2010 © Chiara Modenese <c.modenese@xxxxxxxxx>
-
#include"HertzMindlin.hpp"
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-
YADE_PLUGIN(
(MindlinPhys)
(Ip2_FrictMat_FrictMat_MindlinPhys)
(Law2_ScGeom_MindlinPhys_MindlinDeresiewitz)
(Law2_ScGeom_MindlinPhys_HertzWithLinearShear)
(Law2_ScGeom_MindlinPhys_Mindlin)
+ (MindlinCapillaryPhys)
+ (Ip2_FrictMat_FrictMat_MindlinCapillaryPhys)
);
Real Law2_ScGeom_MindlinPhys_Mindlin::getfrictionDissipation() {return (Real) frictionDissipation;}
@@ -20,12 +20,6 @@
Real Law2_ScGeom_MindlinPhys_Mindlin::getnormDampDissip() {return (Real) normDampDissip;}
Real Law2_ScGeom_MindlinPhys_Mindlin::getshearDampDissip() {return (Real) shearDampDissip;}
-/******************** MindlinPhys *****************************/
-CREATE_LOGGER(MindlinPhys);
-
-MindlinPhys::~MindlinPhys(){}; // destructor
-
-
/******************** Ip2_FrictMat_FrictMat_MindlinPhys *******/
CREATE_LOGGER(Ip2_FrictMat_FrictMat_MindlinPhys);
@@ -36,8 +30,6 @@
FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
-
-
/* from interaction physics */
Real Ea = mat1->young;
Real Eb = mat2->young;
@@ -154,7 +146,6 @@
return adhesionEnergy;
}
-
void Law2_ScGeom_MindlinPhys_MindlinDeresiewitz::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
Body::id_t id1(contact->getId1()), id2(contact->getId2());
ScGeom* geom = static_cast<ScGeom*>(ig.get());
@@ -532,5 +523,74 @@
//phys->prevNormal = scg->normal;
}
+// The following code was moved from Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp
+
+void Ip2_FrictMat_FrictMat_MindlinCapillaryPhys::go( const shared_ptr<Material>& b1 //FrictMat
+ , const shared_ptr<Material>& b2 // FrictMat
+ , const shared_ptr<Interaction>& interaction)
+{
+ if(interaction->phys) return; // no updates of an already existing contact necessary
+
+ shared_ptr<MindlinCapillaryPhys> contactPhysics(new MindlinCapillaryPhys());
+ interaction->phys = contactPhysics;
+
+ FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
+ FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
+
+ /* from interaction physics */
+ Real Ea = mat1->young;
+ Real Eb = mat2->young;
+ Real Va = mat1->poisson;
+ Real Vb = mat2->poisson;
+ Real fa = mat1->frictionAngle;
+ Real fb = mat2->frictionAngle;
+
+ /* from interaction geometry */
+ GenericSpheresContact* scg = YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
+ Real Da = scg->refR1>0 ? scg->refR1 : scg->refR2;
+ Real Db = scg->refR2;
+ //Vector3r normal=scg->normal; //The variable set but not used
+
+ /* calculate stiffness coefficients */
+ Real Ga = Ea/(2*(1+Va));
+ Real Gb = Eb/(2*(1+Vb));
+ Real G = (Ga+Gb)/2; // average of shear modulus
+ Real V = (Va+Vb)/2; // average of poisson's ratio
+ Real E = Ea*Eb/((1.-std::pow(Va,2))*Eb+(1.-std::pow(Vb,2))*Ea); // Young modulus
+ Real R = Da*Db/(Da+Db); // equivalent radius
+ 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 Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force as predicted by DMT theory
+
+ /* pass values calculated from above to MindlinCapillaryPhys */
+ contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
+ //mindlinPhys->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;
+
+ 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_MindlinCapillaryPhys: only one of en, betan can be specified.");
+ if(es && betas) throw std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinCapillaryPhys: only one of es, betas can be specified.");
+
+ // en or es specified, just compute alpha, otherwise alpha remains 0
+ if(en || es){
+ Real logE = log((*en)(mat1->id,mat2->id));
+ 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{
+ 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-09-07 19:40:12 +0000
+++ pkg/dem/HertzMindlin.hpp 2014-07-14 20:08:34 +0000
@@ -18,8 +18,6 @@
#include<yade/pkg/common/NormShearPhys.hpp>
#include<yade/pkg/common/MatchMaker.hpp>
-
-#include <set>
#include <boost/tuple/tuple.hpp>
#include<yade/lib/base/openmp-accu.hpp>
@@ -27,7 +25,7 @@
/******************** MindlinPhys *********************************/
class MindlinPhys: public FrictPhys{
public:
- virtual ~MindlinPhys();
+ virtual ~MindlinPhys() {};
YADE_CLASS_BASE_DOC_ATTRS_CTOR(MindlinPhys,FrictPhys,"Representation of an interaction of the Hertz-Mindlin type.",
((Real,kno,0.0,,"Constant value in the formulation of the normal stiffness"))
((Real,kso,0.0,,"Constant value in the formulation of the tangential stiffness"))
@@ -158,11 +156,53 @@
};
REGISTER_SERIALIZABLE(Law2_ScGeom_MindlinPhys_Mindlin);
-
-
-
-
-
-
-
-
+// The following code was moved from Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp
+
+class MindlinCapillaryPhys : public MindlinPhys
+{
+ public :
+ int currentIndexes [4]; // used for faster interpolation (stores previous positions in tables)
+
+ virtual ~MindlinCapillaryPhys() {};
+
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(MindlinCapillaryPhys,MindlinPhys,"Adds capillary physics to Mindlin's interaction physics.",
+ ((bool,meniscus,false,,"Presence of a meniscus if true"))
+ ((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
+ ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
+ ((Real,vMeniscus,0.,,"Volume of the menicus"))
+ ((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
+ ((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
+ ((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
+ ((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
+ ,/*deprec*/
+ ((Fcap,fCap,"naming convention"))
+ ((CapillaryPressure,capillaryPressure,"naming convention"))
+ ,,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
+ ,
+ );
+ REGISTER_CLASS_INDEX(MindlinCapillaryPhys,MindlinPhys);
+};
+REGISTER_SERIALIZABLE(MindlinCapillaryPhys);
+
+
+class Ip2_FrictMat_FrictMat_MindlinCapillaryPhys : public IPhysFunctor
+{
+ public :
+ virtual void go( const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+
+ FUNCTOR2D(FrictMat,FrictMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_MindlinCapillaryPhys,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
+ ((Real,gamma,0.0,,"Surface energy parameter [J/m^2] per each unit contact surface, to derive DMT formulation from HM"))
+ ((Real,eta,0.0,,"Coefficient to determine the plastic bending moment"))
+ ((Real,krot,0.0,,"Rotational stiffness for moment contact law"))
+ ((Real,ktwist,0.0,,"Torsional stiffness for moment contact law"))
+ ((shared_ptr<MatchMaker>,en,,,"Normal coefficient of restitution $e_n$."))
+ ((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$."))
+ );
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_MindlinCapillaryPhys);
=== removed file 'pkg/dem/InelastCohFrictMat.cpp'
--- pkg/dem/InelastCohFrictMat.cpp 2012-12-14 14:38:44 +0000
+++ pkg/dem/InelastCohFrictMat.cpp 1970-01-01 00:00:00 +0000
@@ -1,7 +0,0 @@
-#include "InelastCohFrictMat.hpp"
-
-InelastCohFrictMat::~InelastCohFrictMat()
-{
-}
-
-YADE_PLUGIN((InelastCohFrictMat));
\ No newline at end of file
=== removed file 'pkg/dem/InelastCohFrictMat.hpp'
--- pkg/dem/InelastCohFrictMat.hpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/InelastCohFrictMat.hpp 1970-01-01 00:00:00 +0000
@@ -1,50 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#pragma once
-
-#include "../common/ElastMat.hpp"
-
-class InelastCohFrictMat : public FrictMat
-{
- public :
- virtual ~InelastCohFrictMat ();
-
-/// Serialization
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictMat,FrictMat,"",
- ((Real,tensionModulus,0.0,,"Tension elasticity modulus"))
- ((Real,compressionModulus,0.0,,"Compresion elasticity modulus"))
- ((Real,shearModulus,0.0,,"shear elasticity modulus"))
- ((Real,alphaKr,2.0,,"Dimensionless coefficient used for the rolling stiffness."))
- ((Real,alphaKtw,2.0,,"Dimensionless coefficient used for the twist stiffness."))
-
- ((Real,nuBending,0.0,,"Bending elastic stress limit"))
- ((Real,nuTwist,0.0,,"Twist elastic stress limit"))
- ((Real,sigmaTension,0.0,,"Tension elastic stress limit"))
- ((Real,sigmaCompression,0.0,,"Compression elastic stress limit"))
- ((Real,shearCohesion,0.0,,"Shear elastic stress limit"))
-
- ((Real,creepTension,0.0,,"Tension/compression creeping coefficient. Usual values between 0 and 1."))
- ((Real,creepBending,0.0,,"Bending creeping coefficient. Usual values between 0 and 1."))
- ((Real,creepTwist,0.0,,"Twist creeping coefficient. Usual values between 0 and 1."))
-
- ((Real,unloadTension,0.0,,"Tension/compression plastic unload coefficient. Usual values between 0 and +infinity."))
- ((Real,unloadBending,0.0,,"Bending plastic unload coefficient. Usual values between 0 and +infinity."))
- ((Real,unloadTwist,0.0,,"Twist plastic unload coefficient. Usual values between 0 and +infinity."))
-
- ((Real,epsilonMaxTension,0.0,,"Maximal plastic strain tension"))
- ((Real,epsilonMaxCompression,0.0,,"Maximal plastic strain compression"))
- ((Real,etaMaxBending,0.0,,"Maximal plastic bending strain"))
- ((Real,etaMaxTwist,0.0,,"Maximal plastic twist strain")),
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(InelastCohFrictMat,FrictMat);
-};
-
-REGISTER_SERIALIZABLE(InelastCohFrictMat);
=== added file 'pkg/dem/InelastCohFrictPM.cpp'
--- pkg/dem/InelastCohFrictPM.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/InelastCohFrictPM.cpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,271 @@
+#include "InelastCohFrictPM.hpp"
+
+YADE_PLUGIN((InelastCohFrictMat)(InelastCohFrictPhys)(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys)(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment));
+
+
+void Ip2_2xInelastCohFrictMat_InelastCohFrictPhys::go(const shared_ptr<Material>& b1 // InelastCohFrictMat
+ , const shared_ptr<Material>& b2 // InelastCohFrictMat
+ , const shared_ptr<Interaction>& interaction)
+{
+
+ InelastCohFrictMat* sdec1 = static_cast<InelastCohFrictMat*>(b1.get());
+ InelastCohFrictMat* sdec2 = static_cast<InelastCohFrictMat*>(b2.get());
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
+
+ //FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
+
+ if (geom) {
+ if (!interaction->phys) {
+ interaction->phys = shared_ptr<InelastCohFrictPhys>(new InelastCohFrictPhys());
+ InelastCohFrictPhys* contactPhysics = YADE_CAST<InelastCohFrictPhys*>(interaction->phys.get());
+ Real pi = 3.14159265;
+ Real r1 = geom->radius1;
+ Real r2 = geom->radius2;
+ Real f1 = sdec1->frictionAngle;
+ Real f2 = sdec2->frictionAngle;
+
+ contactPhysics->tangensOfFrictionAngle = tan(min(f1,f2));
+
+ // harmonic average of modulus
+ contactPhysics->knC = 2.0*sdec1->compressionModulus*r1*sdec2->compressionModulus*r2/(sdec1->compressionModulus*r1+sdec2->compressionModulus*r2);
+ contactPhysics->knT = 2.0*sdec1->tensionModulus*r1*sdec2->tensionModulus*r2/(sdec1->tensionModulus*r1+sdec2->tensionModulus*r2);
+ contactPhysics->ks = 2.0*sdec1->shearModulus*r1*sdec2->shearModulus*r2/(sdec1->shearModulus*r1+sdec2->shearModulus*r2);
+
+ // harmonic average of coeficients for bending and twist coeficients
+ Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
+ Real AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
+
+ contactPhysics->kr = r1*r2*contactPhysics->ks*AlphaKr;
+ contactPhysics->ktw = r1*r2*contactPhysics->ks*AlphaKtw;
+
+ contactPhysics->kTCrp = contactPhysics->knT*min(sdec1->creepTension,sdec2->creepTension);
+ contactPhysics->kRCrp = contactPhysics->kr*min(sdec1->creepBending,sdec2->creepBending);
+ contactPhysics->kTwCrp = contactPhysics->ktw*min(sdec1->creepTwist,sdec2->creepTwist);
+
+ contactPhysics->kRUnld = contactPhysics->kr*min(sdec1->unloadBending,sdec2->unloadBending);
+ contactPhysics->kTUnld = contactPhysics->knT*min(sdec1->unloadTension,sdec2->unloadTension);
+ contactPhysics->kTwUnld = contactPhysics->ktw*min(sdec1->unloadTwist,sdec2->unloadTwist);
+
+ contactPhysics->maxElC = min(sdec1->sigmaCompression,sdec2->sigmaCompression)*pow(min(r2, r1),2);
+ contactPhysics->maxElT = min(sdec1->sigmaTension,sdec2->sigmaTension)*pow(min(r2, r1),2);
+ contactPhysics->maxElB = min(sdec1->nuBending,sdec2->nuBending)*pow(min(r2, r1),3);
+ contactPhysics->maxElTw = min(sdec1->nuTwist,sdec2->nuTwist)*pow(min(r2, r1),3);
+
+ contactPhysics->shearAdhesion = min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(min(r1, r2),2);
+
+ contactPhysics->maxExten = min(sdec1->epsilonMaxTension*r1,sdec2->epsilonMaxTension*r2);
+ contactPhysics->maxContract = min(sdec1->epsilonMaxCompression*r1,sdec2->epsilonMaxCompression*r2);
+
+ contactPhysics->maxBendMom = min(sdec1->etaMaxBending,sdec2->etaMaxBending)*pow(min(r2, r1),3);
+ contactPhysics->maxTwist = 2*pi*min(sdec1->etaMaxTwist,sdec2->etaMaxTwist);
+ }
+ }
+};
+
+
+
+Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy()
+{ //FIXME : this have to be checked and adapted
+ Real normEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
+ if (phys) {
+ normEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
+ }
+ }
+ return normEnergy;
+}
+
+Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy()
+{ //FIXME : this have to be checked and adapted
+ Real shearEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
+ if (phys) {
+ shearEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
+ }
+ }
+ return shearEnergy;
+}
+
+
+void Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
+{
+//FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
+ const int &id1 = contact->getId1();
+ const int &id2 = contact->getId2();
+ const Real& dt = scene->dt;
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*> (ig.get());
+ InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*> (ip.get());
+ if (contact->isFresh(scene)) phys->shearForce = Vector3r::Zero();
+
+ Real un = geom->penetrationDepth-phys->unp;
+ Real Fn;
+
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
+
+
+ if(un<=0){/// tension ///
+ if(-un>phys->maxExten || phys->isBroken){//plastic failure.
+ phys->isBroken=1;
+ phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
+ scene->interactions->requestErase(contact);
+ return;
+ }
+ Fn=phys->knT*un; //elasticity
+ if(-Fn>phys->maxElT || phys->onPlastT){ //so we are on plastic deformation.
+ phys->onPlastT=1;
+ phys->onPlastC=1; //if plasticity is reached on tension, set it to compression too.
+ if(phys->maxCrpRchdT[0]<un){ //unloading/reloading on plastic deformation.
+ Fn = phys->kTUnld*(un-phys->maxCrpRchdT[0])+phys->maxCrpRchdT[1];
+ }
+ else{//loading on plastic deformation : creep.
+ Fn = -phys->maxElT+phys->kTCrp*(un+phys->maxElT/phys->knT);
+ phys->maxCrpRchdT[0]=un; //new maximum is reached.
+ phys->maxCrpRchdT[1]=Fn;
+ }
+ if (Fn>0){ //so the contact just passed the equilibrium state, set new "unp" who stores the plastic equilibrium state.
+ phys->unp=geom->penetrationDepth;
+ phys->maxCrpRchdT[0]=1e20;
+ phys->maxElT=0;
+ }
+ }
+ else{ //elasticity
+ phys->maxCrpRchdT[0]=un;
+ phys->maxCrpRchdT[1]=Fn;
+ }
+ }
+
+ else{/// compression /// similar to tension.
+ if(un>phys->maxContract || phys->isBroken){
+ phys->isBroken=1;
+ phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
+ if(geom->penetrationDepth<=0){ //do not erase the contact while penetrationDepth<0 because it would be recreated at next timestep.
+ scene->interactions->requestErase(contact);
+ }
+ return;
+ }
+ Fn=phys->knC*un;
+ if(Fn>phys->maxElC || phys->onPlastC){
+ phys->onPlastC=1;
+ if(phys->maxCrpRchdC[0]>un){
+ Fn = phys->kTUnld*(un-phys->maxCrpRchdC[0])+phys->maxCrpRchdC[1];
+ }
+ else{
+ Fn = phys->maxElC+phys->kTCrp*(un-phys->maxElC/phys->knC);
+ phys->maxCrpRchdC[0]=un;
+ phys->maxCrpRchdC[1]=Fn;
+ }
+ if (Fn<0){
+ phys->unp=geom->penetrationDepth;
+ phys->maxCrpRchdC[0]=-1e20;
+ phys->maxElC=0;
+ }
+ }
+ else{
+ phys->maxCrpRchdC[0]=un;
+ phys->maxCrpRchdC[1]=Fn;
+ }
+ }
+
+ /// Shear ///
+ Vector3r shearForce = geom->rotate(phys->shearForce);
+ const Vector3r& dus = geom->shearIncrement();
+
+ //Linear elasticity giving "trial" shear force
+ shearForce += phys->ks*dus;
+ Real Fs = shearForce.norm();
+ Real maxFs = phys->shearAdhesion;
+ if (maxFs==0)maxFs = Fn*phys->tangensOfFrictionAngle;
+ maxFs = std::max((Real) 0, maxFs);
+ if (Fs > maxFs) {//Plasticity condition on shear force
+ if (!phys->cohesionBroken) {
+ phys->cohesionBroken=1;
+ phys->shearAdhesion=0;
+ maxFs = max((Real) 0, Fn*phys->tangensOfFrictionAngle);
+ }
+ maxFs = maxFs / Fs;
+ shearForce *= maxFs;
+ }
+
+ //rotational moment are only applied if the cohesion is not broken.
+ /// Twist /// the twist law is driven by twist displacement ("getTwist()").
+ if(!phys->cohesionBroken){
+ Real twist = geom->getTwist() - phys->twp;
+ Real twistM=twist*phys->ktw; //elastic twist moment.
+ bool sgnChanged=0; //whether the twist moment just passed the equilibrium state.
+ if(!contact->isFresh(scene) && phys->moment_twist.dot(twistM*geom->normal)<0)sgnChanged=1;
+ if(std::abs(twist)>phys->maxTwist){
+ phys->cohesionBroken=1;
+ twistM=0;
+ }
+ else{
+ if(std::abs(twistM)>phys->maxElTw || phys->onPlastTw){ //plastic deformation.
+ phys->onPlastTw=1;
+ if(std::abs(phys->maxCrpRchdTw[0])>std::abs(twist)){ //unloading/reloading
+ twistM = phys->kTwUnld*(twist-phys->maxCrpRchdTw[0])+phys->maxCrpRchdTw[1];
+ }
+ else{//creep loading.
+ int sign = twist<0?-1:1;
+ twistM = sign*phys->maxElTw+phys->kTwCrp*(twist-sign*phys->maxElTw/phys->ktw); //creep
+ phys->maxCrpRchdTw[0]=twist; //new maximum reached
+ phys->maxCrpRchdTw[1]=twistM;
+ }
+ if(sgnChanged){
+ phys->maxElTw=0;
+ phys->twp=geom->getTwist();
+ phys->maxCrpRchdTw[0]=0;
+ }
+ }
+ else{ //elasticity
+ phys->maxCrpRchdTw[0]=twist;
+ phys->maxCrpRchdTw[1]=twistM;
+ }
+ }
+ phys->moment_twist = twistM * geom->normal;
+ }
+ else phys->moment_twist=Vector3r(0,0,0);
+
+ /// Bending /// incremental form.
+ if(!phys->cohesionBroken){
+ Vector3r bendM = phys->moment_bending;
+ Vector3r relAngVel = geom->getRelAngVel(de1,de2,dt);
+ Vector3r relRotBend = (relAngVel - geom->normal.dot(relAngVel)*geom->normal)*dt; // relative rotation due to rolling behaviour
+ bendM = geom->rotate(phys->moment_bending); // rotate moment vector (updated)
+ phys->pureCreep=geom->rotate(phys->pureCreep); // pure creep is updated to compute the damage.
+ Vector3r bendM_elast = bendM-phys->kr*relRotBend;
+ if(bendM_elast.norm()>phys->maxElB || phys->onPlastB){ // plastic behavior
+ phys->onPlastB=1;
+ bendM=bendM-phys->kDam*relRotBend; //trial bending
+ if(bendM.norm()<phys->moment_bending.norm()){ // if bending decreased, we are unloading ...
+ bendM = bendM+phys->kDam*relRotBend-phys->kRUnld*relRotBend; // ... so undo bendM and apply unload coefficient.
+ Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend; // trial pure creep.
+ phys->pureCreep = newPureCreep.norm()<phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while unloading, pure creep must decrease.
+ phys->kDam=phys->kr+(phys->kRCrp-phys->kr)*(phys->maxCrpRchdB.norm()-phys->maxElB)/(phys->maxBendMom-phys->maxElB); // compute the damage coefficient.
+ }
+ else{ // bending increased, so we are loading (bendM has to be unchanged).
+ Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend;
+ phys->pureCreep = newPureCreep.norm()>phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while loading, pure creep must increase.
+ if(phys->pureCreep.norm()<bendM.norm()) bendM=phys->pureCreep; // bending moment can't be greather than pure creep.
+ if(phys->pureCreep.norm()>phys->maxCrpRchdB.norm()) phys->maxCrpRchdB=phys->pureCreep; // maxCrpRchdB must follow the maximum of pure creep.
+ if(phys->pureCreep.norm()>phys->maxBendMom){
+ phys->cohesionBroken=1;
+ bendM=bendM_elast=Vector3r(0,0,0);
+ }
+ }
+ phys->moment_bending=bendM;
+ }
+ else{//elasticity
+ phys->pureCreep=phys->moment_bending=phys->maxCrpRchdB=bendM_elast;
+ phys->kDam=phys->kRCrp;
+ }
+ }
+ phys->shearForce=shearForce;
+ phys->normalForce=-Fn*geom->normal;
+ applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
+ scene->forces.addTorque(id1,-phys->moment_bending-phys->moment_twist);
+ scene->forces.addTorque(id2,phys->moment_bending+phys->moment_twist);
+}
=== added file 'pkg/dem/InelastCohFrictPM.hpp'
--- pkg/dem/InelastCohFrictPM.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/InelastCohFrictPM.hpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,149 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
+#pragma once
+
+#include <yade/pkg/common/ElastMat.hpp>
+#include <yade/pkg/dem/ScGeom.hpp>
+#include <yade/pkg/dem/FrictPhys.hpp>
+#include "CohesiveFrictionalContactLaw.hpp"
+
+class InelastCohFrictMat : public FrictMat
+{
+ public :
+ virtual ~InelastCohFrictMat () {};
+
+/// Serialization
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictMat,FrictMat,"",
+ ((Real,tensionModulus,0.0,,"Tension elasticity modulus"))
+ ((Real,compressionModulus,0.0,,"Compresion elasticity modulus"))
+ ((Real,shearModulus,0.0,,"shear elasticity modulus"))
+ ((Real,alphaKr,2.0,,"Dimensionless coefficient used for the rolling stiffness."))
+ ((Real,alphaKtw,2.0,,"Dimensionless coefficient used for the twist stiffness."))
+
+ ((Real,nuBending,0.0,,"Bending elastic stress limit"))
+ ((Real,nuTwist,0.0,,"Twist elastic stress limit"))
+ ((Real,sigmaTension,0.0,,"Tension elastic stress limit"))
+ ((Real,sigmaCompression,0.0,,"Compression elastic stress limit"))
+ ((Real,shearCohesion,0.0,,"Shear elastic stress limit"))
+
+ ((Real,creepTension,0.0,,"Tension/compression creeping coefficient. Usual values between 0 and 1."))
+ ((Real,creepBending,0.0,,"Bending creeping coefficient. Usual values between 0 and 1."))
+ ((Real,creepTwist,0.0,,"Twist creeping coefficient. Usual values between 0 and 1."))
+
+ ((Real,unloadTension,0.0,,"Tension/compression plastic unload coefficient. Usual values between 0 and +infinity."))
+ ((Real,unloadBending,0.0,,"Bending plastic unload coefficient. Usual values between 0 and +infinity."))
+ ((Real,unloadTwist,0.0,,"Twist plastic unload coefficient. Usual values between 0 and +infinity."))
+
+ ((Real,epsilonMaxTension,0.0,,"Maximal plastic strain tension"))
+ ((Real,epsilonMaxCompression,0.0,,"Maximal plastic strain compression"))
+ ((Real,etaMaxBending,0.0,,"Maximal plastic bending strain"))
+ ((Real,etaMaxTwist,0.0,,"Maximal plastic twist strain")),
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(InelastCohFrictMat,FrictMat);
+};
+
+REGISTER_SERIALIZABLE(InelastCohFrictMat);
+
+class InelastCohFrictPhys : public FrictPhys
+{
+ public :
+ virtual ~InelastCohFrictPhys() {};
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictPhys,FrictPhys,"",
+ ((bool,cohesionBroken,false,,"is cohesion active? will be set false when a fragile contact is broken"))
+
+ ((Real,knT,0,,"tension stiffness"))
+ ((Real,knC,0,,"compression stiffness"))
+ ((Real,ktw,0,,"twist shear stiffness"))
+ ((Real,ks,0,,"shear stiffness"))
+ ((Real,kr,0,,"bending stiffness"))
+
+ ((Real,maxElB,0.0,,"Maximum bending elastic moment."))
+ ((Real,maxElTw,0.0,,"Maximum twist elastic moment."))
+ ((Real,maxElT,0.0,,"Maximum tension elastic force."))
+ ((Real,maxElC,0.0,,"Maximum compression elastic force."))
+ ((Real,shearAdhesion,0,,"Maximum elastic shear force (cohesion)."))
+
+ ((Real,kTCrp,0.0,,"Tension/compression creep stiffness"))
+ ((Real,kRCrp,0.0,,"Bending creep stiffness"))
+ ((Real,kTwCrp,0.0,,"Twist creep stiffness"))
+
+ ((Real,kTUnld,0.0,,"Tension/compression plastic unload stiffness"))
+ ((Real,kRUnld,0.0,,"Bending plastic unload stiffness"))
+ ((Real,kTwUnld,0.0,,"Twist plastic unload stiffness"))
+
+ ((Real,maxExten,0.0,,"Plastic failure extension (stretching)."))
+ ((Real,maxContract,0.0,,"Plastic failure contraction (shrinkage)."))
+ ((Real,maxBendMom,0.0,,"Plastic failure bending moment."))
+ ((Real,maxTwist,0.0,,"Plastic failure twist angle"))
+
+ ((bool,isBroken,false,,"true if compression plastic fracture achieved"))
+
+ ((Real,unp,0,,"plastic normal penetration depth describing the equilibrium state."))
+ ((Real,twp,0,,"plastic twist penetration depth describing the equilibrium state."))
+
+ ((bool,onPlastB,false,Attr::readonly,"true if plasticity achieved on bending"))
+ ((bool,onPlastTw,false,Attr::readonly,"true if plasticity achieved on twisting"))
+ ((bool,onPlastT,false,Attr::readonly,"true if plasticity achieved on traction"))
+ ((bool,onPlastC,false,Attr::readonly,"true if plasticity achieved on compression"))
+
+ ((Vector2r,maxCrpRchdT,Vector2r(0,0),Attr::readonly,"maximal extension reached on plastic deformation. maxCrpRchdT[0] stores un and maxCrpRchdT[1] stores Fn."))
+ ((Vector2r,maxCrpRchdC,Vector2r(0,0),Attr::readonly,"maximal compression reached on plastic deformation. maxCrpRchdC[0] stores un and maxCrpRchdC[1] stores Fn."))
+ ((Vector2r,maxCrpRchdTw,Vector2r(0,0),Attr::readonly,"maximal twist reached on plastic deformation. maxCrpRchdTw[0] stores twist angle and maxCrpRchdTw[1] stores twist moment."))
+ ((Vector3r,maxCrpRchdB,Vector3r(0,0,0),Attr::readonly,"maximal bending moment reached on plastic deformation."))
+
+ ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::readonly),"Twist moment"))
+ ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::readonly),"Bending moment"))
+ ((Vector3r,pureCreep,Vector3r(0,0,0),(Attr::readonly),"Pure creep curve, used for comparison in calculation."))
+ ((Real,kDam,0,(Attr::readonly),"Damage coefficient on bending, computed from maximum bending moment reached and pure creep behaviour. Its values will vary between :yref:`InelastCohFrictPhys::kr` and :yref:`InelastCohFrictPhys::kRCrp` ."))
+ // internal attributes
+ ,
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(InelastCohFrictPhys,FrictPhys);
+
+};
+
+REGISTER_SERIALIZABLE(InelastCohFrictPhys);
+
+class Ip2_2xInelastCohFrictMat_InelastCohFrictPhys : public IPhysFunctor
+{
+ public :
+ virtual void go(const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+ int cohesionDefinitionIteration;
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys,IPhysFunctor,
+ "Generates cohesive-frictional interactions with moments. Used in the contact law :yref:`Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment`.",
+ ,
+ cohesionDefinitionIteration = -1;
+ );
+ FUNCTOR2D(InelastCohFrictMat,InelastCohFrictMat);
+};
+
+REGISTER_SERIALIZABLE(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys);
+
+
+class Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment: public LawFunctor{
+ public:
+ Real normElastEnergy();
+ Real shearElastEnergy();
+ virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment,LawFunctor,"This law is currently under developpement. Final version and documentation will come before the end of 2014.",
+ ,,
+ .def("normElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
+ .def("shearElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
+ );
+ FUNCTOR2D(ScGeom6D,InelastCohFrictPhys);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment);
=== removed file 'pkg/dem/InelastCohFrictPhys.cpp'
--- pkg/dem/InelastCohFrictPhys.cpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/InelastCohFrictPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,13 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#include "InelastCohFrictPhys.hpp"
-InelastCohFrictPhys::~InelastCohFrictPhys()
-{
-}
-YADE_PLUGIN((InelastCohFrictPhys));
\ No newline at end of file
=== removed file 'pkg/dem/InelastCohFrictPhys.hpp'
--- pkg/dem/InelastCohFrictPhys.hpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/InelastCohFrictPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,74 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#pragma once
-
-#include "FrictPhys.hpp"
-
-
-class InelastCohFrictPhys : public FrictPhys
-{
- public :
- virtual ~InelastCohFrictPhys();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictPhys,FrictPhys,"",
- ((bool,cohesionBroken,false,,"is cohesion active? will be set false when a fragile contact is broken"))
-
- ((Real,knT,0,,"tension stiffness"))
- ((Real,knC,0,,"compression stiffness"))
- ((Real,ktw,0,,"twist shear stiffness"))
- ((Real,ks,0,,"shear stiffness"))
- ((Real,kr,0,,"bending stiffness"))
-
- ((Real,maxElB,0.0,,"Maximum bending elastic moment."))
- ((Real,maxElTw,0.0,,"Maximum twist elastic moment."))
- ((Real,maxElT,0.0,,"Maximum tension elastic force."))
- ((Real,maxElC,0.0,,"Maximum compression elastic force."))
- ((Real,shearAdhesion,0,,"Maximum elastic shear force (cohesion)."))
-
- ((Real,kTCrp,0.0,,"Tension/compression creep stiffness"))
- ((Real,kRCrp,0.0,,"Bending creep stiffness"))
- ((Real,kTwCrp,0.0,,"Twist creep stiffness"))
-
- ((Real,kTUnld,0.0,,"Tension/compression plastic unload stiffness"))
- ((Real,kRUnld,0.0,,"Bending plastic unload stiffness"))
- ((Real,kTwUnld,0.0,,"Twist plastic unload stiffness"))
-
- ((Real,maxExten,0.0,,"Plastic failure extension (stretching)."))
- ((Real,maxContract,0.0,,"Plastic failure contraction (shrinkage)."))
- ((Real,maxBendMom,0.0,,"Plastic failure bending moment."))
- ((Real,maxTwist,0.0,,"Plastic failure twist angle"))
-
- ((bool,isBroken,false,,"true if compression plastic fracture achieved"))
-
- ((Real,unp,0,,"plastic normal penetration depth describing the equilibrium state."))
- ((Real,twp,0,,"plastic twist penetration depth describing the equilibrium state."))
-
- ((bool,onPlastB,false,Attr::readonly,"true if plasticity achieved on bending"))
- ((bool,onPlastTw,false,Attr::readonly,"true if plasticity achieved on twisting"))
- ((bool,onPlastT,false,Attr::readonly,"true if plasticity achieved on traction"))
- ((bool,onPlastC,false,Attr::readonly,"true if plasticity achieved on compression"))
-
- ((Vector2r,maxCrpRchdT,Vector2r(0,0),Attr::readonly,"maximal extension reached on plastic deformation. maxCrpRchdT[0] stores un and maxCrpRchdT[1] stores Fn."))
- ((Vector2r,maxCrpRchdC,Vector2r(0,0),Attr::readonly,"maximal compression reached on plastic deformation. maxCrpRchdC[0] stores un and maxCrpRchdC[1] stores Fn."))
- ((Vector2r,maxCrpRchdTw,Vector2r(0,0),Attr::readonly,"maximal twist reached on plastic deformation. maxCrpRchdTw[0] stores twist angle and maxCrpRchdTw[1] stores twist moment."))
- ((Vector3r,maxCrpRchdB,Vector3r(0,0,0),Attr::readonly,"maximal bending moment reached on plastic deformation."))
-
- ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::readonly),"Twist moment"))
- ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::readonly),"Bending moment"))
- ((Vector3r,pureCreep,Vector3r(0,0,0),(Attr::readonly),"Pure creep curve, used for comparison in calculation."))
- ((Real,kDam,0,(Attr::readonly),"Damage coefficient on bending, computed from maximum bending moment reached and pure creep behaviour. Its values will vary between :yref:`InelastCohFrictPhys::kr` and :yref:`InelastCohFrictPhys::kRCrp` ."))
- // internal attributes
- ,
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(InelastCohFrictPhys,FrictPhys);
-
-};
-
-REGISTER_SERIALIZABLE(InelastCohFrictPhys);
=== removed file 'pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp'
--- pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,72 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#include "Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-
-
-
-void Ip2_2xInelastCohFrictMat_InelastCohFrictPhys::go(const shared_ptr<Material>& b1 // InelastCohFrictMat
- , const shared_ptr<Material>& b2 // InelastCohFrictMat
- , const shared_ptr<Interaction>& interaction)
-{
-
- InelastCohFrictMat* sdec1 = static_cast<InelastCohFrictMat*>(b1.get());
- InelastCohFrictMat* sdec2 = static_cast<InelastCohFrictMat*>(b2.get());
- ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
-
- //FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
-
- if (geom) {
- if (!interaction->phys) {
- interaction->phys = shared_ptr<InelastCohFrictPhys>(new InelastCohFrictPhys());
- InelastCohFrictPhys* contactPhysics = YADE_CAST<InelastCohFrictPhys*>(interaction->phys.get());
- Real pi = 3.14159265;
- Real r1 = geom->radius1;
- Real r2 = geom->radius2;
- Real f1 = sdec1->frictionAngle;
- Real f2 = sdec2->frictionAngle;
-
- contactPhysics->tangensOfFrictionAngle = tan(min(f1,f2));
-
- // harmonic average of modulus
- contactPhysics->knC = 2.0*sdec1->compressionModulus*r1*sdec2->compressionModulus*r2/(sdec1->compressionModulus*r1+sdec2->compressionModulus*r2);
- contactPhysics->knT = 2.0*sdec1->tensionModulus*r1*sdec2->tensionModulus*r2/(sdec1->tensionModulus*r1+sdec2->tensionModulus*r2);
- contactPhysics->ks = 2.0*sdec1->shearModulus*r1*sdec2->shearModulus*r2/(sdec1->shearModulus*r1+sdec2->shearModulus*r2);
-
- // harmonic average of coeficients for bending and twist coeficients
- Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
- Real AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
-
- contactPhysics->kr = r1*r2*contactPhysics->ks*AlphaKr;
- contactPhysics->ktw = r1*r2*contactPhysics->ks*AlphaKtw;
-
- contactPhysics->kTCrp = contactPhysics->knT*min(sdec1->creepTension,sdec2->creepTension);
- contactPhysics->kRCrp = contactPhysics->kr*min(sdec1->creepBending,sdec2->creepBending);
- contactPhysics->kTwCrp = contactPhysics->ktw*min(sdec1->creepTwist,sdec2->creepTwist);
-
- contactPhysics->kRUnld = contactPhysics->kr*min(sdec1->unloadBending,sdec2->unloadBending);
- contactPhysics->kTUnld = contactPhysics->knT*min(sdec1->unloadTension,sdec2->unloadTension);
- contactPhysics->kTwUnld = contactPhysics->ktw*min(sdec1->unloadTwist,sdec2->unloadTwist);
-
- contactPhysics->maxElC = min(sdec1->sigmaCompression,sdec2->sigmaCompression)*pow(min(r2, r1),2);
- contactPhysics->maxElT = min(sdec1->sigmaTension,sdec2->sigmaTension)*pow(min(r2, r1),2);
- contactPhysics->maxElB = min(sdec1->nuBending,sdec2->nuBending)*pow(min(r2, r1),3);
- contactPhysics->maxElTw = min(sdec1->nuTwist,sdec2->nuTwist)*pow(min(r2, r1),3);
-
- contactPhysics->shearAdhesion = min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(min(r1, r2),2);
-
- contactPhysics->maxExten = min(sdec1->epsilonMaxTension*r1,sdec2->epsilonMaxTension*r2);
- contactPhysics->maxContract = min(sdec1->epsilonMaxCompression*r1,sdec2->epsilonMaxCompression*r2);
-
- contactPhysics->maxBendMom = min(sdec1->etaMaxBending,sdec2->etaMaxBending)*pow(min(r2, r1),3);
- contactPhysics->maxTwist = 2*pi*min(sdec1->etaMaxTwist,sdec2->etaMaxTwist);
- }
- }
-};
-YADE_PLUGIN((Ip2_2xInelastCohFrictMat_InelastCohFrictPhys));
\ No newline at end of file
=== removed file 'pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp'
--- pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,32 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#pragma once
-#include "Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp"
-#include "InelastCohFrictMat.hpp"
-#include "InelastCohFrictPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-
-
-class Ip2_2xInelastCohFrictMat_InelastCohFrictPhys : public IPhysFunctor
-{
- public :
- virtual void go( const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
- int cohesionDefinitionIteration;
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys,IPhysFunctor,
- "Generates cohesive-frictional interactions with moments. Used in the contact law :yref:`Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment`.",
- ,
- cohesionDefinitionIteration = -1;
- );
- FUNCTOR2D(InelastCohFrictMat,InelastCohFrictMat);
-};
-
-REGISTER_SERIALIZABLE(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys);
\ No newline at end of file
=== removed file 'pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp'
--- pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 2011-08-26 09:49:25 +0000
+++ pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,69 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/NormalInelasticityPhys.hpp>
-#include<yade/pkg/dem/NormalInelasticMat.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-
-
-
-void Ip2_2xNormalInelasticMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // NormalInelasticMat
- , const shared_ptr<Material>& b2 // NormalInelasticMat
- , const shared_ptr<Interaction>& interaction)
-{
- NormalInelasticMat* sdec1 = static_cast<NormalInelasticMat*>(b1.get());
- NormalInelasticMat* sdec2 = static_cast<NormalInelasticMat*>(b2.get());
- ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
-
-
- if(geom) // so it is ScGeom - NON PERMANENT LINK
- {
- if(!interaction->phys)
- {
-//std::cerr << " isNew, id1: " << interaction->getId1() << " id2: " << interaction->getId2() << "\n";
- interaction->phys = shared_ptr<NormalInelasticityPhys>(new NormalInelasticityPhys());
- NormalInelasticityPhys* contactPhysics = YADE_CAST<NormalInelasticityPhys*>(interaction->phys.get());
-
- Real Ea = sdec1->young;
- Real Eb = sdec2->young;
- Real Va = sdec1->poisson;
- Real Vb = sdec2->poisson;
- Real Ra = geom->radius1;
- Real Rb = geom->radius2;
- Real fa = sdec1->frictionAngle;
- Real fb = sdec2->frictionAngle;
-
- Real Kn = 2.0*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);//harmonic average of two stiffnesses
-
- Real Ks = 2.0*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
-
- // Jean-Patrick Plassiard, Noura Belheine, Frederic Victor Donze, "A Spherical Discrete Element Model: calibration procedure and incremental response", DOI: 10.1007/s10035-009-0130-x
-
- Real Kr = betaR*std::pow((Ra+Rb)/2.0,2)*Ks;
-
- contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
- contactPhysics->forMaxMoment = 1.0*(Ra+Rb)/2.0; // 1.0 corresponding to ethaR which I don't know exactly where to define as a parameter...
-
- // Lot of suppress here around (>) r2276.
- contactPhysics->knLower = Kn;
- contactPhysics->kn = Kn;
- contactPhysics->ks = Ks;
- contactPhysics->kr = Kr;
- }
-
- }
-
-};
-YADE_PLUGIN((Ip2_2xNormalInelasticMat_NormalInelasticityPhys));
-
-
-
=== removed file 'pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp'
--- pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp 2011-05-04 10:19:03 +0000
+++ pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,38 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/common/Dispatching.hpp>
-#include<yade/pkg/dem/NormalInelasticMat.hpp>
-
-/*! \brief The RelationShips for using Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity
-
-In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed.
-WARNING : as in the others Relationships most of the attributes are computed only once : when the interaction is "new"
- */
-
-class Ip2_2xNormalInelasticMat_NormalInelasticityPhys : public IPhysFunctor
-{
- public :
-
- virtual void go( const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
-
- FUNCTOR2D(NormalInelasticMat,NormalInelasticMat);
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_2xNormalInelasticMat_NormalInelasticityPhys,
- IPhysFunctor,
- "The RelationShips for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`\n\n In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
- ((Real,betaR,0.12,,"Parameter for computing the torque-stifness : T-stifness = betaR * Rmoy^2"))
- );
-};
-
-REGISTER_SERIALIZABLE(Ip2_2xNormalInelasticMat_NormalInelasticityPhys);
-
-
=== removed file 'pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp 2014-05-26 16:49:23 +0000
+++ pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,89 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/CohFrictPhys.hpp>
-#include<yade/pkg/dem/CohFrictMat.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-void Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::go(const shared_ptr<Material>& b1 // CohFrictMat
- , const shared_ptr<Material>& b2 // CohFrictMat
- , const shared_ptr<Interaction>& interaction)
-{
- CohFrictMat* sdec1 = static_cast<CohFrictMat*>(b1.get());
- CohFrictMat* sdec2 = static_cast<CohFrictMat*>(b2.get());
- ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
-
- //Create cohesive interractions only once
- if (setCohesionNow && cohesionDefinitionIteration==-1) cohesionDefinitionIteration=scene->iter;
- if (setCohesionNow && cohesionDefinitionIteration!=-1 && cohesionDefinitionIteration!=scene->iter) {
- cohesionDefinitionIteration = -1;
- setCohesionNow = 0;}
-
- if (geom) {
- if (!interaction->phys) {
- interaction->phys = shared_ptr<CohFrictPhys>(new CohFrictPhys());
- CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());
- Real Ea = sdec1->young;
- Real Eb = sdec2->young;
- Real Va = sdec1->poisson;
- Real Vb = sdec2->poisson;
- Real Da = geom->radius1;
- Real Db = geom->radius2;
- Real fa = sdec1->frictionAngle;
- Real fb = sdec2->frictionAngle;
- Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
-
- // harmonic average of alphas parameters
- Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
- Real AlphaKtw;
- if (sdec1->alphaKtw && sdec2->alphaKtw) AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
- else AlphaKtw=0;
-
- Real Ks;
- if (Va && Vb) Ks = 2.0*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
- else Ks=0;
-
- contactPhysics->kr = Da*Db*Ks*AlphaKr;
- contactPhysics->ktw = Da*Db*Ks*AlphaKtw;
- contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
-
- if ((setCohesionOnNewContacts || setCohesionNow) && sdec1->isCohesive && sdec2->isCohesive)
- {
- contactPhysics->cohesionBroken = false;
- contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(Db, Da),2);
- contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(Db, Da),2);
- geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
- }
- contactPhysics->kn = Kn;
- contactPhysics->ks = Ks;
-
- contactPhysics->maxRollPl = min(sdec1->etaRoll*Da,sdec2->etaRoll*Db);
- contactPhysics->maxTwistPl = min(sdec1->etaTwist*Da,sdec2->etaTwist*Db);
- contactPhysics->momentRotationLaw=(sdec1->momentRotationLaw && sdec2->momentRotationLaw);
- }
- else {// !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
- CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());
- if ((setCohesionNow && sdec1->isCohesive && sdec2->isCohesive) || contactPhysics->initCohesion)
- {
- contactPhysics->cohesionBroken = false;
- contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
- contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
-
- geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
- contactPhysics->initCohesion=false;
- }
- }
- }
-};
-YADE_PLUGIN((Ip2_CohFrictMat_CohFrictMat_CohFrictPhys));
-
-
-
=== removed file 'pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp'
--- pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp 2014-05-21 16:01:55 +0000
+++ pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,34 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
-* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/common/Dispatching.hpp>
-#include<yade/pkg/dem/CohFrictMat.hpp>
-
-class Ip2_CohFrictMat_CohFrictMat_CohFrictPhys : public IPhysFunctor
-{
- public :
- virtual void go( const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
- int cohesionDefinitionIteration;
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_CohFrictMat_CohFrictMat_CohFrictPhys,IPhysFunctor,
- "Generates cohesive-frictional interactions with moments, used in the contact law :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`. The normal/shear stiffness and friction definitions are the same as in :yref:`Ip2_FrictMat_FrictMat_FrictPhys`, check the documentation there for details.\n\nAdhesions related to the normal and the shear components are calculated form :yref:`CohFrictMat::normalCohesion` ($C_n$) and :yref:`CohFrictMat::shearlCohesion` ($C_s$). For particles of size $R_1$,$R_2$ the adhesion will be $a_i=C_i min(R_1,R_2)^2$, $i=n\\,s$.\n\nTwist and rolling stiffnesses are proportional to the shear stiffness through dimensionless factors alphaKtw and alphaKr, such that the rotational stiffnesses are defined by $k_s \\alpha_i R_1 R_2$, $i=tw\\,r$",
- ((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts in current time-step. The flag is turned false automatically, so that assignment is done in the current timestep only."))
- ((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts. If false, only existing contacts can be cohesive (also see :yref:`Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::setCohesionNow`), and new contacts are only frictional."))
- ,
- cohesionDefinitionIteration = -1;
- );
- FUNCTOR2D(CohFrictMat,CohFrictMat);
-};
-
-REGISTER_SERIALIZABLE(Ip2_CohFrictMat_CohFrictMat_CohFrictPhys);
-
-
=== modified file 'pkg/dem/Ip2_ElastMat.cpp'
--- pkg/dem/Ip2_ElastMat.cpp 2013-09-23 17:39:59 +0000
+++ pkg/dem/Ip2_ElastMat.cpp 2014-07-14 20:08:34 +0000
@@ -1,7 +1,6 @@
#include "Ip2_ElastMat.hpp"
#include <yade/pkg/common/NormShearPhys.hpp>
-#include <yade/pkg/common/NormShearPhys.hpp>
#include <yade/pkg/dem/DemXDofGeom.hpp>
YADE_PLUGIN((Ip2_ElastMat_ElastMat_NormPhys)(Ip2_ElastMat_ElastMat_NormShearPhys));
=== removed file 'pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp 2010-12-31 14:35:21 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,56 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno CHAREYRE *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"Ip2_FrictMat_FrictMat_CapillaryPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/CapillaryPhys.hpp>
-#include<yade/pkg/dem/FrictPhys.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-#include<yade/pkg/common/Dispatching.hpp>
-
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-YADE_PLUGIN((Ip2_FrictMat_FrictMat_CapillaryPhys));
-
-void Ip2_FrictMat_FrictMat_CapillaryPhys::go( const shared_ptr<Material>& b1 //FrictMat
- , const shared_ptr<Material>& b2 // FrictMat
- , const shared_ptr<Interaction>& interaction)
-{
- ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
- if(geom)
- {
- if(!interaction->phys)
- {
- const shared_ptr<FrictMat>& sdec1 = YADE_PTR_CAST<FrictMat>(b1);
- const shared_ptr<FrictMat>& sdec2 = YADE_PTR_CAST<FrictMat>(b2);
-
- if (!interaction->phys) interaction->phys = shared_ptr<CapillaryPhys>(new CapillaryPhys());
- const shared_ptr<CapillaryPhys>& contactPhysics = YADE_PTR_CAST<CapillaryPhys>(interaction->phys);
-
- Real Ea = sdec1->young;
- Real Eb = sdec2->young;
- Real Va = sdec1->poisson;
- Real Vb = sdec2->poisson;
- Real Da = geom->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size)
- Real Db = geom->radius2; // FIXME - as above
- Real fa = sdec1->frictionAngle;
- Real fb = sdec2->frictionAngle;
- Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
- Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
-
- contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
- contactPhysics->kn = Kn;
- contactPhysics->ks = Ks;
- }
- }
-};
-
-
-
-
=== removed file 'pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.hpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.hpp 2012-03-27 22:07:35 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_CapillaryPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,29 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno CHAREYRE *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/common/Dispatching.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-
-
-class Ip2_FrictMat_FrictMat_CapillaryPhys : public IPhysFunctor
-{
- public :
- virtual void go( const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
-
- FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_FrictMat_FrictMat_CapillaryPhys,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",,;);
-
-};
-REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_CapillaryPhys);
-
-
-
=== removed file 'pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2014-01-23 10:25:43 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,86 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno CHAREYRE *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"Ip2_FrictMat_FrictMat_FrictPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/DemXDofGeom.hpp>
-#include<yade/pkg/dem/FrictPhys.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-#include<yade/pkg/dem/ViscoelasticPM.hpp>
-#include <cassert>
-
-
-
-void Ip2_FrictMat_FrictMat_FrictPhys::go( const shared_ptr<Material>& b1
- , const shared_ptr<Material>& b2
- , const shared_ptr<Interaction>& interaction)
-{
- if(interaction->phys) return;
-
- const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
- const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
-
- Real Ra,Rb;//Vector3r normal;
- assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode
- GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
- Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
- Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
-
- interaction->phys = shared_ptr<FrictPhys>(new FrictPhys());
- const shared_ptr<FrictPhys>& contactPhysics = YADE_PTR_CAST<FrictPhys>(interaction->phys);
- Real Ea = mat1->young;
- Real Eb = mat2->young;
- Real Va = mat1->poisson;
- Real Vb = mat2->poisson;
-
- //harmonic average of the two stiffnesses when (2*Ri*Ei) is the stiffness of a contact point on sphere "i"
- Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
- //same for shear stiffness
- Real Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);
-
- Real frictionAngle = (!frictAngle) ? std::min(mat1->frictionAngle,mat2->frictionAngle) : (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle);
- contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
- contactPhysics->kn = Kn;
- contactPhysics->ks = Ks;
-};
-YADE_PLUGIN((Ip2_FrictMat_FrictMat_FrictPhys));
-
-void Ip2_FrictMat_FrictMat_ViscoFrictPhys::go( const shared_ptr<Material>& b1
- , const shared_ptr<Material>& b2
- , const shared_ptr<Interaction>& interaction)
-{
- if(interaction->phys) return;
- const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
- const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
- interaction->phys = shared_ptr<ViscoFrictPhys>(new ViscoFrictPhys());
- const shared_ptr<ViscoFrictPhys>& contactPhysics = YADE_PTR_CAST<ViscoFrictPhys>(interaction->phys);
- Real Ea = mat1->young;
- Real Eb = mat2->young;
- Real Va = mat1->poisson;
- Real Vb = mat2->poisson;
-
- Real Ra,Rb;//Vector3r normal;
- assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode
- GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
- Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
- Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
-
- //harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
- Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
- //same for shear stiffness
- Real Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);
-
- Real frictionAngle = (!frictAngle) ? std::min(mat1->frictionAngle,mat2->frictionAngle) : (*frictAngle)(mat1->id,mat2->id,mat1->frictionAngle,mat2->frictionAngle);
- contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
- contactPhysics->kn = Kn;
- contactPhysics->ks = Ks;
-};
-YADE_PLUGIN((Ip2_FrictMat_FrictMat_ViscoFrictPhys));
-
=== removed file 'pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp 2014-01-23 10:25:43 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,38 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno CHAREYRE *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/common/Dispatching.hpp>
-#include<yade/pkg/common/MatchMaker.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-
-class Ip2_FrictMat_FrictMat_FrictPhys: public IPhysFunctor{
- public:
- virtual void go(const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
- FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_FrictPhys,IPhysFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under point load is defined here as $1/(E.D)$, with $E$ the stiffness of the sphere and $D$ its diameter. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. $1/(E_1.D_1)+1/(E_2.D_2)$ in the general case, or $2/(E.D)$ in the special case of equal sizes and equal stiffness. Note that summing compliances corresponds to an harmonic average of stiffnesss (as in e.g. [Scholtes2009a]_), which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor:\n\n $k_n = \\frac{E_1D_1*E_2D_2}{E_1D_1+E_2D_2}=\\frac{k_1*k_2}{k_1+k_2}$, with $k_i=E_iD_i$.\n\n The shear stiffness ks of one sphere is defined via the material parameter :yref:`ElastMat::poisson`, as ks=poisson*kn, and the resulting shear stiffness of the interaction will be also an harmonic average. In the case of a contact between a :yref:`ViscElMat` and a :yref:`FrictMat`, be sure to set :yref:`FrictMat::young` and :yref:`FrictMat::poisson`, otherwise the default value will be used.",
- ((shared_ptr<MatchMaker>,frictAngle,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's friction angle. If ``None``, minimum value is used."))
- );
-};
-REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_FrictPhys);
-
-class Ip2_FrictMat_FrictMat_ViscoFrictPhys: public Ip2_FrictMat_FrictMat_FrictPhys{
- public:
- virtual void go(const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
- FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_ViscoFrictPhys,Ip2_FrictMat_FrictMat_FrictPhys,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under symetric point loads is defined here as 1/(E.r), with E the stiffness of the sphere and r its radius, and corresponds to a compliance 1/(2.E.r)=1/(E.D) from each contact point. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. 1/(E.D1)+1/(E.D2) in the general case, or 1/(E.r) in the special case of equal sizes. Note that summing compliances corresponds to an harmonic average of stiffnesss, which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor.\n\nThe shear stiffness ks of one sphere is defined via the material parameter :yref:`ElastMat::poisson`, as ks=poisson*kn, and the resulting shear stiffness of the interaction will be also an harmonic average.",
- );
-};
-REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_ViscoFrictPhys);
-
-
=== removed file 'pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp 2012-02-16 16:05:15 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,95 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno CHAREYRE *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include"Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/HertzMindlin.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-#include<yade/pkg/common/Dispatching.hpp>
-
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-YADE_PLUGIN(
- (MindlinCapillaryPhys)
- (Ip2_FrictMat_FrictMat_MindlinCapillaryPhys)
-);
-
-MindlinCapillaryPhys::~MindlinCapillaryPhys(){};// destructor
-
-void Ip2_FrictMat_FrictMat_MindlinCapillaryPhys::go( const shared_ptr<Material>& b1 //FrictMat
- , const shared_ptr<Material>& b2 // FrictMat
- , const shared_ptr<Interaction>& interaction)
-{
- if(interaction->phys) return; // no updates of an already existing contact necessary
-
- shared_ptr<MindlinCapillaryPhys> contactPhysics(new MindlinCapillaryPhys());
- interaction->phys = contactPhysics;
-
- FrictMat* mat1 = YADE_CAST<FrictMat*>(b1.get());
- FrictMat* mat2 = YADE_CAST<FrictMat*>(b2.get());
-
- /* from interaction physics */
- Real Ea = mat1->young;
- Real Eb = mat2->young;
- Real Va = mat1->poisson;
- Real Vb = mat2->poisson;
- Real fa = mat1->frictionAngle;
- Real fb = mat2->frictionAngle;
-
- /* from interaction geometry */
- GenericSpheresContact* scg = YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
- Real Da = scg->refR1>0 ? scg->refR1 : scg->refR2;
- Real Db = scg->refR2;
- //Vector3r normal=scg->normal; //The variable set but not used
-
- /* calculate stiffness coefficients */
- Real Ga = Ea/(2*(1+Va));
- Real Gb = Eb/(2*(1+Vb));
- Real G = (Ga+Gb)/2; // average of shear modulus
- Real V = (Va+Vb)/2; // average of poisson's ratio
- Real E = Ea*Eb/((1.-std::pow(Va,2))*Eb+(1.-std::pow(Vb,2))*Ea); // Young modulus
- Real R = Da*Db/(Da+Db); // equivalent radius
- 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 Adhesion = 4.*Mathr::PI*R*gamma; // calculate adhesion force as predicted by DMT theory
-
- /* pass values calculated from above to MindlinCapillaryPhys */
- contactPhysics->tangensOfFrictionAngle = std::tan(frictionAngle);
- //mindlinPhys->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;
-
- 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_MindlinCapillaryPhys: only one of en, betan can be specified.");
- if(es && betas) throw std::invalid_argument("Ip2_FrictMat_FrictMat_MindlinCapillaryPhys: only one of es, betas can be specified.");
-
- // en or es specified, just compute alpha, otherwise alpha remains 0
- if(en || es){
- Real logE = log((*en)(mat1->id,mat2->id));
- 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{
- contactPhysics->betan=betan ? (*betan)(mat1->id,mat2->id) : 0;
- contactPhysics->betas=betas ? (*betas)(mat1->id,mat2->id) : contactPhysics->betan;
- }
-};
-
-
-
-
=== removed file 'pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp 2013-05-30 20:17:45 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,65 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2007 by Bruno CHAREYRE *
-* bruno.chareyre@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/common/Dispatching.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-#include<yade/pkg/dem/HertzMindlin.hpp>
-
-class MindlinCapillaryPhys : public MindlinPhys
-{
- public :
- int currentIndexes [4]; // used for faster interpolation (stores previous positions in tables)
-
- virtual ~MindlinCapillaryPhys();
-
- YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(MindlinCapillaryPhys,MindlinPhys,"Adds capillary physics to Mindlin's interaction physics.",
- ((bool,meniscus,false,,"Presence of a meniscus if true"))
- ((bool,isBroken,false,,"If true, capillary force is zero and liquid bridge is inactive."))
- ((Real,capillaryPressure,0.,,"Value of the capillary pressure Uc defines as Ugas-Uliquid"))
- ((Real,vMeniscus,0.,,"Volume of the menicus"))
- ((Real,Delta1,0.,,"Defines the surface area wetted by the meniscus on the smallest grains of radius R1 (R1<R2)"))
- ((Real,Delta2,0.,,"Defines the surface area wetted by the meniscus on the biggest grains of radius R2 (R1<R2)"))
- ((Vector3r,fCap,Vector3r::Zero(),,"Capillary Force produces by the presence of the meniscus"))
- ((short int,fusionNumber,0.,,"Indicates the number of meniscii that overlap with this one"))
- ,/*deprec*/
- ((Fcap,fCap,"naming convention"))
- ((CapillaryPressure,capillaryPressure,"naming convention"))
- ,,createIndex();currentIndexes[0]=currentIndexes[1]=currentIndexes[2]=currentIndexes[3]=0;
- ,
- );
- REGISTER_CLASS_INDEX(MindlinCapillaryPhys,MindlinPhys);
-};
-REGISTER_SERIALIZABLE(MindlinCapillaryPhys);
-
-
-class Ip2_FrictMat_FrictMat_MindlinCapillaryPhys : public IPhysFunctor
-{
- public :
- virtual void go( const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
-
- FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictMat_MindlinCapillaryPhys,IPhysFunctor, "RelationShips to use with Law2_ScGeom_CapillaryPhys_Capillarity\n\n In these RelationShips all the interaction attributes are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
- ((Real,gamma,0.0,,"Surface energy parameter [J/m^2] per each unit contact surface, to derive DMT formulation from HM"))
- ((Real,eta,0.0,,"Coefficient to determine the plastic bending moment"))
- ((Real,krot,0.0,,"Rotational stiffness for moment contact law"))
- ((Real,ktwist,0.0,,"Torsional stiffness for moment contact law"))
- ((shared_ptr<MatchMaker>,en,,,"Normal coefficient of restitution $e_n$."))
- ((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$."))
- );
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_MindlinCapillaryPhys);
-
-
-
=== removed file 'pkg/dem/KinemCNDEngine.cpp'
--- pkg/dem/KinemCNDEngine.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/KinemCNDEngine.cpp 1970-01-01 00:00:00 +0000
@@ -1,53 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#include<yade/pkg/dem/KinemCNDEngine.hpp>
-
-
-YADE_PLUGIN((KinemCNDEngine));
-
-void KinemCNDEngine::action()
-{
- KinemSimpleShearBox::getBoxes_Dt();
-
- if( ((shearSpeed > 0) && (gamma<=gammalim)) || ((shearSpeed < 0) /*&& (gamma>=gammalim)*/ ) )
- {
- if(temoinfin!=0)
- temoinfin=0;
- letMove(shearSpeed * dt,0);
- gamma+=shearSpeed * dt;
- }
- else
- {
- stopMovement();
- if(temoinfin==0)
- {
- Omega::instance().saveSimulation(Key + "endShear.xml");
- temoinfin=1;
- }
- }
-
- for(unsigned int j=0;j<gamma_save.size();j++)
- {
- if ( ( ( (shearSpeed>0)&&(gamma > gamma_save[j]) ) || ((shearSpeed<0)&&(gamma < gamma_save[j])) ) && (temoin_save[j]==0) )
- {
- stopMovement(); // reset of all the speeds before the save
- Omega::instance().saveSimulation(Key+"_"+boost::lexical_cast<string> (floor(gamma*1000)) + "mmsheared.xml");
- temoin_save[j]=1;
- }
- }
-
-}
-
-
-
-
-
-
-
=== removed file 'pkg/dem/KinemCNDEngine.hpp'
--- pkg/dem/KinemCNDEngine.hpp 2011-04-12 22:06:51 +0000
+++ pkg/dem/KinemCNDEngine.hpp 1970-01-01 00:00:00 +0000
@@ -1,38 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/dem/KinemSimpleShearBox.hpp>
-
-
-class KinemCNDEngine : public KinemSimpleShearBox
-{
- private :
- int temoinfin;
-
-
-
- public :
- void action();
-
- protected :
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCNDEngine,KinemSimpleShearBox,
- "To apply a Constant Normal Displacement (CND) shear for a parallelogram box\n\n \tThis engine, designed for simulations implying a simple shear box (:yref:`SimpleShear` Preprocessor or scripts/simpleShear.py), allows one to perform a constant normal displacement shear, by translating horizontally the upper plate, while the lateral ones rotate so that they always keep contact with the lower and upper walls.",
- ((Real,shearSpeed,0.0,,"the speed at which the shear is performed : speed of the upper plate [m/s]"))
- ((Real,gammalim,0.0,,"the value of the tangential displacement at wich the displacement is stopped [m]"))
- ((Real,gamma,0.0,,"the current value of the tangential displacement"))
- ((std::vector<Real>,gamma_save,,,"vector with the values of gamma at which a save of the simulation is performed [m]")),
- temoinfin=0;
- );
-};
-
-REGISTER_SERIALIZABLE(KinemCNDEngine);
-
-
=== removed file 'pkg/dem/KinemCNLEngine.cpp'
--- pkg/dem/KinemCNLEngine.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/KinemCNLEngine.cpp 1970-01-01 00:00:00 +0000
@@ -1,61 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include<yade/pkg/dem/KinemCNLEngine.hpp>
-
-
-YADE_PLUGIN((KinemCNLEngine));
-
-void KinemCNLEngine::action()
-{
- if(LOG) cout << "debut applyCondi du CNCEngine !!" << endl;
- KinemSimpleShearBox::getBoxes_Dt();
-
- if(LOG) cout << "gamma = " << boost::lexical_cast<string>(gamma) << " et gammalim = " << boost::lexical_cast<string>(gammalim) << endl;
- if(gamma<=gammalim)
- {
- if(LOG) cout << "Je suis bien dans la partie gamma < gammalim" << endl;
- if(temoin==0)
-
- {
- if(LOG) cout << "Je veux maintenir la Force a f0 = : " << f0 << endl;
- temoin=1;
- }
- computeDY(0.0);
-
- letMove(shearSpeed*dt,deltaH);
- gamma+=shearSpeed * dt;
-
- }
- else if (temoin<2)
- {
- stopMovement(); // INDISPENSABLE !
- it_stop=scene->iter;
- cout << "Shear stopped : gammaLim reached at it "<< it_stop << endl;
- temoin=2;
- }
- else if (temoin==2 && (scene->iter==(it_stop+5000)) )
- {
- Omega::instance().saveSimulation(Key + "endShear" +boost::lexical_cast<string> ( scene->iter ) + ".xml");
- Omega::instance().pause();
- }
-
- for(unsigned int j=0;j<gamma_save.size();j++)
- {
- if ((gamma > gamma_save[j]) && (temoin_save[j]==0))
- {
- stopMovement(); // reset of all the speeds before the save
- Omega::instance().saveSimulation(Key+"_"+boost::lexical_cast<string> (floor(gamma*1000)) +"_" +boost::lexical_cast<string> (floor(gamma*10000)-10*floor(gamma*1000))+ "mmsheared.xml");
- temoin_save[j]=1;
- }
- }
-
-}
-
-
-
=== removed file 'pkg/dem/KinemCNLEngine.hpp'
--- pkg/dem/KinemCNLEngine.hpp 2011-04-12 22:06:51 +0000
+++ pkg/dem/KinemCNLEngine.hpp 1970-01-01 00:00:00 +0000
@@ -1,43 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/dem/KinemSimpleShearBox.hpp>
-
-
-
-
-class KinemCNLEngine : public KinemSimpleShearBox
-{
- private :
-
- int temoin,// utile pour savoir ou on en est
- it_stop
- ;
-
- public :
- void action()
- ;
-
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCNLEngine,KinemSimpleShearBox,
- "To apply a constant normal stress shear (i.e. Constant Normal Load : CNL) for a parallelogram box (simple shear box : :yref:`SimpleShear` Preprocessor or scripts/simpleShear.py)\n\nThis engine allows one to translate horizontally the upper plate while the lateral ones rotate so that they always keep contact with the lower and upper walls.\n\nIn fact the upper plate can move not only horizontally but also vertically, so that the normal stress acting on it remains constant (this constant value is not chosen by the user but is the one that exists at the beginning of the simulation)\n\nThe right vertical displacements which will be allowed are computed from the rigidity Kn of the sample over the wall (so to cancel a deltaSigma, a normal dplt deltaSigma*S/(Kn) is set)\n\nThe movement is moreover controlled by the user via a *shearSpeed* which will be the speed of the upper wall, and by a maximum value of horizontal displacement *gammalim*, after which the shear stops.\n\n.. note::\n\tNot only the positions of walls are updated but also their speeds, which is all but useless considering the fact that in the contact laws these velocities of bodies are used to compute values of tangential relative displacements.\n\n.. warning::\n\tBecause of this last point, if you want to use later saves of simulations executed with this Engine, but without that stopMovement was executed, your boxes will keep their speeds => you will have to cancel them 'by hand' in the .xml.\n",
- ((Real,shearSpeed,0.0,,"the speed at wich the shearing is performed : speed of the upper plate [m/s]"))
- ((Real,gammalim,0.0,,"the value of tangential displacement (of upper plate) at wich the shearing is stopped [m]"))
- ((Real,gamma,0.0,,"current value of tangential displacement [m]"))
- ((std::vector<Real>,gamma_save,,,"vector with the values of gamma at which a save of the simulation is performed [m]")),
-// ((Real,coeff_dech,1.0,,"in the case of the use of 'Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity' for ex, where kn(unload)#kn(load). The engine cares to find the value at the first run BROKEN actually")),
- temoin=0;
- it_stop=0;
- );
-};
-
-REGISTER_SERIALIZABLE(KinemCNLEngine);
-
-
=== removed file 'pkg/dem/KinemCNSEngine.cpp'
--- pkg/dem/KinemCNSEngine.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/KinemCNSEngine.cpp 1970-01-01 00:00:00 +0000
@@ -1,51 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-
-#include<yade/pkg/dem/KinemCNSEngine.hpp>
-
-
-YADE_PLUGIN((KinemCNSEngine));
-
-
-void KinemCNSEngine::action()
-{
- if(LOG) cerr << "debut applyCondi !!" << endl;
- KinemSimpleShearBox::getBoxes_Dt();
-
- if(gamma<=gammalim)
- {
- computeDY(KnC);
- letMove(shearSpeed * dt,deltaH);
- gamma+=shearSpeed * dt;
- if(temoin==0)
- {
- temoin=1;
- }
- }
- else if (temoin<2)
- {
- stopMovement(); // INDISPENSABLE !
- it_stop=scene->iter;
- cout << "Cisaillement arrete : gammaLim atteint a l'iteration "<< it_stop << endl;
- temoin=2;
- }
- else if (temoin==2 && (scene->iter==(it_stop+5000)) )
- {
- Omega::instance().saveSimulation(Key + "finCis" +boost::lexical_cast<string> (scene->iter ) + ".xml");
- Omega::instance().pause();
- }
-
-}
-
-
-
-
-
-
=== removed file 'pkg/dem/KinemCNSEngine.hpp'
--- pkg/dem/KinemCNSEngine.hpp 2011-04-12 22:06:51 +0000
+++ pkg/dem/KinemCNSEngine.hpp 1970-01-01 00:00:00 +0000
@@ -1,41 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/dem/KinemSimpleShearBox.hpp>
-
-
-
-class KinemCNSEngine : public KinemSimpleShearBox
-{
- private :
- int temoin,it_stop
- ;
-
- public :
- void action()
- ;
-
-
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCNSEngine,KinemSimpleShearBox,
- "To apply a Constant Normal Stifness (CNS) shear for a parallelogram box (simple shear)\n\nThis engine, useable in simulations implying one deformable parallelepipedic box, allows one to translate horizontally the upper plate while the lateral ones rotate so that they always keep contact with the lower and upper walls. The upper plate can move not only horizontally but also vertically, so that the normal rigidity defined by DeltaF(upper plate)/DeltaU(upper plate) = constant (= *KnC* defined by the user).\n\nThe movement is moreover controlled by the user via a *shearSpeed* which is the horizontal speed of the upper wall, and by a maximum value of horizontal displacement *gammalim* (of the upper plate), after which the shear stops.\n\n.. note::\n\t not only the positions of walls are updated but also their speeds, which is all but useless considering the fact that in the contact laws these velocities of bodies are used to compute values of tangential relative displacements.\n\n.. warning::\n\tBut, because of this last point, if you want to use later saves of simulations executed with this Engine, but without that stopMovement was executed, your boxes will keep their speeds => you will have to cancel them by hand in the .xml",
- ((Real,shearSpeed,0.0,,"the speed at wich the shearing is performed : speed of the upper plate [m/s]"))
- ((Real,gammalim,0.0,,"the value of tangential displacement (of upper plate) at wich the shearing is stopped [m]"))
- ((Real,gamma,0.0,,"current value of tangential displacement [m]"))
- ((Real,KnC,10.0e6,,"the normal rigidity chosen by the user [MPa/mm] - the conversion in Pa/m will be made")),
- temoin=0;
- it_stop=0;
- );
-
-};
-
-REGISTER_SERIALIZABLE(KinemCNSEngine);
-
-
=== removed file 'pkg/dem/KinemCTDEngine.cpp'
--- pkg/dem/KinemCTDEngine.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/KinemCTDEngine.cpp 1970-01-01 00:00:00 +0000
@@ -1,73 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include<yade/pkg/dem/KinemCTDEngine.hpp>
-
-
-YADE_PLUGIN((KinemCTDEngine));
-
-KinemCTDEngine::~KinemCTDEngine()
-{
-}
-
-void KinemCTDEngine::action()
-{
- KinemSimpleShearBox::getBoxes_Dt();
-
-
- scene->forces.sync();
- Real current_NormalForce=(scene->forces.getForce(id_topbox)).y();
- KinemSimpleShearBox::computeScontact();
- current_sigma=current_NormalForce/(1000.0*Scontact); // so we have the current value of sigma, in kPa
-
- if( ((compSpeed > 0) && (current_sigma < targetSigma)) || ((compSpeed < 0) && (current_sigma > targetSigma)) )
- {
- if(temoin!=0)
- {
-// cout << "j'ai ici un temoin #0 visiblement. En effet temoin =" <<boost::lexical_cast<string>(temoin) << endl;
- temoin=0;
-// cout << "Maintenant (toujours dans le if temoin!=0), temoin =" <<boost::lexical_cast<string>(temoin) << endl;
- }
-
- letMove(0.0,-compSpeed*dt);
- }
- else if (temoin==0)
- {
- stopMovement();
-// cout << "Mouvement stoppe, temoin = " << boost::lexical_cast<string>(temoin) << endl;
-// cout << " Dans le if, temoin =" << boost::lexical_cast<string>(temoin) << endl;
- string f;
- if (compSpeed > 0)
- f="Sigmax_";
- else
- f="Sigmin_";
-
- Omega::instance().saveSimulation(Key + f +boost::lexical_cast<string> (floor(targetSigma)) + "kPaReached.xml");
- temoin=1;
-// cout << " Fin du if, temoin =" << boost::lexical_cast<string>(temoin) << endl << endl;
- }
-
-
- for(unsigned int j=0;j<sigma_save.size();j++)
- {
- if( ( ( (compSpeed>0)&&(current_sigma > sigma_save[j]) ) || ((compSpeed<0)&&(current_sigma < sigma_save[j])) ) && (temoin_save[j]==0))
- {
- stopMovement();
- Omega::instance().saveSimulation(Key + "SigInt_" +boost::lexical_cast<string> (floor(current_sigma)) + "kPareached.xml");
- temoin_save[j]=1;
- }
- }
-// cout << "Fin de ApplyCondi, temoin = " << boost::lexical_cast<string>(temoin) << endl;
-
-}
-
-
-
-
-
-
=== removed file 'pkg/dem/KinemCTDEngine.hpp'
--- pkg/dem/KinemCTDEngine.hpp 2013-05-29 09:48:51 +0000
+++ pkg/dem/KinemCTDEngine.hpp 1970-01-01 00:00:00 +0000
@@ -1,41 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jerome Duriez *
-* jerome.duriez@xxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/dem/KinemSimpleShearBox.hpp>
-
-
-class KinemCTDEngine : public KinemSimpleShearBox
-{
- private :
- Real current_sigma // Computed in kPa
- ;
-
- int temoin;
-
- public :
- virtual ~KinemCTDEngine();
-
- void action();
-// ;
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCTDEngine,KinemSimpleShearBox,
- "To compress a simple shear sample by moving the upper box in a vertical way only, so that the tangential displacement (defined by the horizontal gap between the upper and lower boxes) remains constant (thus, the CTD = Constant Tangential Displacement).\n \t The lateral boxes move also to keep always contact. All that until this box is submitted to a given stress (=*targetSigma*). Moreover saves are executed at each value of stresses stored in the vector *sigma_save*, and at *targetSigma*",
- ((Real,compSpeed,0.0,,"(vertical) speed of the upper box : >0 for real compression, <0 for unloading [$m/s$]"))
- ((std::vector<Real>,sigma_save,,,"vector with the values of sigma at which a save of the simulation should be performed [$kPa$]"))
- ((Real,targetSigma,0.0,,"the value of sigma at which the compression should stop [$kPa$]")),
- temoin=0;
- )
-
-
-};
-
-REGISTER_SERIALIZABLE(KinemCTDEngine);
-
-
=== added file 'pkg/dem/KinemC__Engine.cpp'
--- pkg/dem/KinemC__Engine.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/KinemC__Engine.cpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,169 @@
+/*************************************************************************
+* Copyright (C) 2008 by Jerome Duriez *
+* jerome.duriez@xxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#include<yade/pkg/dem/KinemC__Engine.hpp>
+
+YADE_PLUGIN((KinemCTDEngine)(KinemCNDEngine)(KinemCNLEngine)(KinemCNSEngine));
+
+void KinemCTDEngine::action()
+{
+ KinemSimpleShearBox::getBoxes_Dt();
+
+
+ scene->forces.sync();
+ Real current_NormalForce=(scene->forces.getForce(id_topbox)).y();
+ KinemSimpleShearBox::computeScontact();
+ current_sigma=current_NormalForce/(1000.0*Scontact); // so we have the current value of sigma, in kPa
+
+ if( ((compSpeed > 0) && (current_sigma < targetSigma)) || ((compSpeed < 0) && (current_sigma > targetSigma)) )
+ {
+ if(temoin!=0)
+ {
+// cout << "j'ai ici un temoin #0 visiblement. En effet temoin =" <<boost::lexical_cast<string>(temoin) << endl;
+ temoin=0;
+// cout << "Maintenant (toujours dans le if temoin!=0), temoin =" <<boost::lexical_cast<string>(temoin) << endl;
+ }
+
+ letMove(0.0,-compSpeed*dt);
+ }
+ else if (temoin==0)
+ {
+ stopMovement();
+// cout << "Mouvement stoppe, temoin = " << boost::lexical_cast<string>(temoin) << endl;
+// cout << " Dans le if, temoin =" << boost::lexical_cast<string>(temoin) << endl;
+ string f;
+ if (compSpeed > 0)
+ f="Sigmax_";
+ else
+ f="Sigmin_";
+
+ Omega::instance().saveSimulation(Key + f +boost::lexical_cast<string> (floor(targetSigma)) + "kPaReached.xml");
+ temoin=1;
+// cout << " Fin du if, temoin =" << boost::lexical_cast<string>(temoin) << endl << endl;
+ }
+
+
+ for(unsigned int j=0;j<sigma_save.size();j++)
+ {
+ if( ( ( (compSpeed>0)&&(current_sigma > sigma_save[j]) ) || ((compSpeed<0)&&(current_sigma < sigma_save[j])) ) && (temoin_save[j]==0))
+ {
+ stopMovement();
+ Omega::instance().saveSimulation(Key + "SigInt_" +boost::lexical_cast<string> (floor(current_sigma)) + "kPareached.xml");
+ temoin_save[j]=1;
+ }
+ }
+// cout << "Fin de ApplyCondi, temoin = " << boost::lexical_cast<string>(temoin) << endl;
+
+}
+
+void KinemCNDEngine::action()
+{
+ KinemSimpleShearBox::getBoxes_Dt();
+
+ if( ((shearSpeed > 0) && (gamma<=gammalim)) || ((shearSpeed < 0) /*&& (gamma>=gammalim)*/ ) )
+ {
+ if(temoinfin!=0)
+ temoinfin=0;
+ letMove(shearSpeed * dt,0);
+ gamma+=shearSpeed * dt;
+ }
+ else
+ {
+ stopMovement();
+ if(temoinfin==0)
+ {
+ Omega::instance().saveSimulation(Key + "endShear.xml");
+ temoinfin=1;
+ }
+ }
+
+ for(unsigned int j=0;j<gamma_save.size();j++)
+ {
+ if ( ( ( (shearSpeed>0)&&(gamma > gamma_save[j]) ) || ((shearSpeed<0)&&(gamma < gamma_save[j])) ) && (temoin_save[j]==0) )
+ {
+ stopMovement(); // reset of all the speeds before the save
+ Omega::instance().saveSimulation(Key+"_"+boost::lexical_cast<string> (floor(gamma*1000)) + "mmsheared.xml");
+ temoin_save[j]=1;
+ }
+ }
+
+}
+
+void KinemCNLEngine::action()
+{
+ if(LOG) cout << "debut applyCondi du CNCEngine !!" << endl;
+ KinemSimpleShearBox::getBoxes_Dt();
+
+ if(LOG) cout << "gamma = " << boost::lexical_cast<string>(gamma) << " et gammalim = " << boost::lexical_cast<string>(gammalim) << endl;
+ if(gamma<=gammalim)
+ {
+ if(LOG) cout << "Je suis bien dans la partie gamma < gammalim" << endl;
+ if(temoin==0)
+
+ {
+ if(LOG) cout << "Je veux maintenir la Force a f0 = : " << f0 << endl;
+ temoin=1;
+ }
+ computeDY(0.0);
+
+ letMove(shearSpeed*dt,deltaH);
+ gamma+=shearSpeed * dt;
+
+ }
+ else if (temoin<2)
+ {
+ stopMovement(); // INDISPENSABLE !
+ it_stop=scene->iter;
+ cout << "Shear stopped : gammaLim reached at it "<< it_stop << endl;
+ temoin=2;
+ }
+ else if (temoin==2 && (scene->iter==(it_stop+5000)) )
+ {
+ Omega::instance().saveSimulation(Key + "endShear" +boost::lexical_cast<string> ( scene->iter ) + ".xml");
+ Omega::instance().pause();
+ }
+
+ for(unsigned int j=0;j<gamma_save.size();j++)
+ {
+ if ((gamma > gamma_save[j]) && (temoin_save[j]==0))
+ {
+ stopMovement(); // reset of all the speeds before the save
+ Omega::instance().saveSimulation(Key+"_"+boost::lexical_cast<string> (floor(gamma*1000)) +"_" +boost::lexical_cast<string> (floor(gamma*10000)-10*floor(gamma*1000))+ "mmsheared.xml");
+ temoin_save[j]=1;
+ }
+ }
+}
+
+void KinemCNSEngine::action()
+{
+ if(LOG) cerr << "debut applyCondi !!" << endl;
+ KinemSimpleShearBox::getBoxes_Dt();
+
+ if(gamma<=gammalim)
+ {
+ computeDY(KnC);
+ letMove(shearSpeed * dt,deltaH);
+ gamma+=shearSpeed * dt;
+ if(temoin==0)
+ {
+ temoin=1;
+ }
+ }
+ else if (temoin<2)
+ {
+ stopMovement(); // INDISPENSABLE !
+ it_stop=scene->iter;
+ cout << "Cisaillement arrete : gammaLim atteint a l'iteration "<< it_stop << endl;
+ temoin=2;
+ }
+ else if (temoin==2 && (scene->iter==(it_stop+5000)) )
+ {
+ Omega::instance().saveSimulation(Key + "finCis" +boost::lexical_cast<string> (scene->iter ) + ".xml");
+ Omega::instance().pause();
+ }
+}
=== added file 'pkg/dem/KinemC__Engine.hpp'
--- pkg/dem/KinemC__Engine.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/KinemC__Engine.hpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,108 @@
+/*************************************************************************
+* Copyright (C) 2008 by Jerome Duriez *
+* jerome.duriez@xxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#pragma once
+
+#include<yade/pkg/dem/KinemSimpleShearBox.hpp>
+
+class KinemCTDEngine : public KinemSimpleShearBox
+{
+ private :
+ Real current_sigma // Computed in kPa
+ ;
+
+ int temoin;
+
+ public :
+ virtual ~KinemCTDEngine() {};
+
+ void action();
+// ;
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCTDEngine,KinemSimpleShearBox,
+ "To compress a simple shear sample by moving the upper box in a vertical way only, so that the tangential displacement (defined by the horizontal gap between the upper and lower boxes) remains constant (thus, the CTD = Constant Tangential Displacement).\n \t The lateral boxes move also to keep always contact. All that until this box is submitted to a given stress (=*targetSigma*). Moreover saves are executed at each value of stresses stored in the vector *sigma_save*, and at *targetSigma*",
+ ((Real,compSpeed,0.0,,"(vertical) speed of the upper box : >0 for real compression, <0 for unloading [$m/s$]"))
+ ((std::vector<Real>,sigma_save,,,"vector with the values of sigma at which a save of the simulation should be performed [$kPa$]"))
+ ((Real,targetSigma,0.0,,"the value of sigma at which the compression should stop [$kPa$]")),
+ temoin=0;
+ )
+
+
+};
+
+REGISTER_SERIALIZABLE(KinemCTDEngine);
+
+class KinemCNDEngine : public KinemSimpleShearBox
+{
+ private :
+ int temoinfin;
+ public :
+ void action();
+
+ protected :
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCNDEngine,KinemSimpleShearBox,
+ "To apply a Constant Normal Displacement (CND) shear for a parallelogram box\n\n \tThis engine, designed for simulations implying a simple shear box (:yref:`SimpleShear` Preprocessor or scripts/simpleShear.py), allows one to perform a constant normal displacement shear, by translating horizontally the upper plate, while the lateral ones rotate so that they always keep contact with the lower and upper walls.",
+ ((Real,shearSpeed,0.0,,"the speed at which the shear is performed : speed of the upper plate [m/s]"))
+ ((Real,gammalim,0.0,,"the value of the tangential displacement at wich the displacement is stopped [m]"))
+ ((Real,gamma,0.0,,"the current value of the tangential displacement"))
+ ((std::vector<Real>,gamma_save,,,"vector with the values of gamma at which a save of the simulation is performed [m]")),
+ temoinfin=0;
+ );
+};
+
+REGISTER_SERIALIZABLE(KinemCNDEngine);
+
+class KinemCNLEngine : public KinemSimpleShearBox
+{
+ private :
+
+ int temoin,// utile pour savoir ou on en est
+ it_stop
+ ;
+
+ public :
+ void action()
+ ;
+
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCNLEngine,KinemSimpleShearBox,
+ "To apply a constant normal stress shear (i.e. Constant Normal Load : CNL) for a parallelogram box (simple shear box : :yref:`SimpleShear` Preprocessor or scripts/simpleShear.py)\n\nThis engine allows one to translate horizontally the upper plate while the lateral ones rotate so that they always keep contact with the lower and upper walls.\n\nIn fact the upper plate can move not only horizontally but also vertically, so that the normal stress acting on it remains constant (this constant value is not chosen by the user but is the one that exists at the beginning of the simulation)\n\nThe right vertical displacements which will be allowed are computed from the rigidity Kn of the sample over the wall (so to cancel a deltaSigma, a normal dplt deltaSigma*S/(Kn) is set)\n\nThe movement is moreover controlled by the user via a *shearSpeed* which will be the speed of the upper wall, and by a maximum value of horizontal displacement *gammalim*, after which the shear stops.\n\n.. note::\n\tNot only the positions of walls are updated but also their speeds, which is all but useless considering the fact that in the contact laws these velocities of bodies are used to compute values of tangential relative displacements.\n\n.. warning::\n\tBecause of this last point, if you want to use later saves of simulations executed with this Engine, but without that stopMovement was executed, your boxes will keep their speeds => you will have to cancel them 'by hand' in the .xml.\n",
+ ((Real,shearSpeed,0.0,,"the speed at wich the shearing is performed : speed of the upper plate [m/s]"))
+ ((Real,gammalim,0.0,,"the value of tangential displacement (of upper plate) at wich the shearing is stopped [m]"))
+ ((Real,gamma,0.0,,"current value of tangential displacement [m]"))
+ ((std::vector<Real>,gamma_save,,,"vector with the values of gamma at which a save of the simulation is performed [m]")),
+// ((Real,coeff_dech,1.0,,"in the case of the use of 'Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity' for ex, where kn(unload)#kn(load). The engine cares to find the value at the first run BROKEN actually")),
+ temoin=0;
+ it_stop=0;
+ );
+};
+
+REGISTER_SERIALIZABLE(KinemCNLEngine);
+
+class KinemCNSEngine : public KinemSimpleShearBox
+{
+ private :
+ int temoin,it_stop
+ ;
+ public :
+ void action()
+ ;
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(KinemCNSEngine,KinemSimpleShearBox,
+ "To apply a Constant Normal Stifness (CNS) shear for a parallelogram box (simple shear)\n\nThis engine, useable in simulations implying one deformable parallelepipedic box, allows one to translate horizontally the upper plate while the lateral ones rotate so that they always keep contact with the lower and upper walls. The upper plate can move not only horizontally but also vertically, so that the normal rigidity defined by DeltaF(upper plate)/DeltaU(upper plate) = constant (= *KnC* defined by the user).\n\nThe movement is moreover controlled by the user via a *shearSpeed* which is the horizontal speed of the upper wall, and by a maximum value of horizontal displacement *gammalim* (of the upper plate), after which the shear stops.\n\n.. note::\n\t not only the positions of walls are updated but also their speeds, which is all but useless considering the fact that in the contact laws these velocities of bodies are used to compute values of tangential relative displacements.\n\n.. warning::\n\tBut, because of this last point, if you want to use later saves of simulations executed with this Engine, but without that stopMovement was executed, your boxes will keep their speeds => you will have to cancel them by hand in the .xml",
+ ((Real,shearSpeed,0.0,,"the speed at wich the shearing is performed : speed of the upper plate [m/s]"))
+ ((Real,gammalim,0.0,,"the value of tangential displacement (of upper plate) at wich the shearing is stopped [m]"))
+ ((Real,gamma,0.0,,"current value of tangential displacement [m]"))
+ ((Real,KnC,10.0e6,,"the normal rigidity chosen by the user [MPa/mm] - the conversion in Pa/m will be made")),
+ temoin=0;
+ it_stop=0;
+ );
+
+};
+
+REGISTER_SERIALIZABLE(KinemCNSEngine);
=== removed file 'pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp'
--- pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp 1970-01-01 00:00:00 +0000
@@ -1,226 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-
-
-#include "Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp"
-
-Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy()
-{ //FIXME : this have to be checked and adapted
- Real normEnergy=0;
- FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
- if(!I->isReal()) continue;
- InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
- if (phys) {
- normEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
- }
- }
- return normEnergy;
-}
-Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy()
-{ //FIXME : this have to be checked and adapted
- Real shearEnergy=0;
- FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
- if(!I->isReal()) continue;
- InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*>(I->phys.get());
- if (phys) {
- shearEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
- }
- }
- return shearEnergy;
-}
-
-
-
-void Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
-{
-//FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
- const int &id1 = contact->getId1();
- const int &id2 = contact->getId2();
- const Real& dt = scene->dt;
- ScGeom6D* geom = YADE_CAST<ScGeom6D*> (ig.get());
- InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*> (ip.get());
- if (contact->isFresh(scene)) phys->shearForce = Vector3r::Zero();
-
- Real un = geom->penetrationDepth-phys->unp;
- Real Fn;
-
- State* de1 = Body::byId(id1,scene)->state.get();
- State* de2 = Body::byId(id2,scene)->state.get();
-
-
- if(un<=0){/// tension ///
- if(-un>phys->maxExten || phys->isBroken){//plastic failure.
- phys->isBroken=1;
- phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
- scene->interactions->requestErase(contact);
- return;
- }
- Fn=phys->knT*un; //elasticity
- if(-Fn>phys->maxElT || phys->onPlastT){ //so we are on plastic deformation.
- phys->onPlastT=1;
- phys->onPlastC=1; //if plasticity is reached on tension, set it to compression too.
- if(phys->maxCrpRchdT[0]<un){ //unloading/reloading on plastic deformation.
- Fn = phys->kTUnld*(un-phys->maxCrpRchdT[0])+phys->maxCrpRchdT[1];
- }
- else{//loading on plastic deformation : creep.
- Fn = -phys->maxElT+phys->kTCrp*(un+phys->maxElT/phys->knT);
- phys->maxCrpRchdT[0]=un; //new maximum is reached.
- phys->maxCrpRchdT[1]=Fn;
- }
- if (Fn>0){ //so the contact just passed the equilibrium state, set new "unp" who stores the plastic equilibrium state.
- phys->unp=geom->penetrationDepth;
- phys->maxCrpRchdT[0]=1e20;
- phys->maxElT=0;
- }
- }
- else{ //elasticity
- phys->maxCrpRchdT[0]=un;
- phys->maxCrpRchdT[1]=Fn;
- }
- }
-
- else{/// compression /// similar to tension.
- if(un>phys->maxContract || phys->isBroken){
- phys->isBroken=1;
- phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
- if(geom->penetrationDepth<=0){ //do not erase the contact while penetrationDepth<0 because it would be recreated at next timestep.
- scene->interactions->requestErase(contact);
- }
- return;
- }
- Fn=phys->knC*un;
- if(Fn>phys->maxElC || phys->onPlastC){
- phys->onPlastC=1;
- if(phys->maxCrpRchdC[0]>un){
- Fn = phys->kTUnld*(un-phys->maxCrpRchdC[0])+phys->maxCrpRchdC[1];
- }
- else{
- Fn = phys->maxElC+phys->kTCrp*(un-phys->maxElC/phys->knC);
- phys->maxCrpRchdC[0]=un;
- phys->maxCrpRchdC[1]=Fn;
- }
- if (Fn<0){
- phys->unp=geom->penetrationDepth;
- phys->maxCrpRchdC[0]=-1e20;
- phys->maxElC=0;
- }
- }
- else{
- phys->maxCrpRchdC[0]=un;
- phys->maxCrpRchdC[1]=Fn;
- }
- }
-
- /// Shear ///
- Vector3r shearForce = geom->rotate(phys->shearForce);
- const Vector3r& dus = geom->shearIncrement();
-
- //Linear elasticity giving "trial" shear force
- shearForce += phys->ks*dus;
- Real Fs = shearForce.norm();
- Real maxFs = phys->shearAdhesion;
- if (maxFs==0)maxFs = Fn*phys->tangensOfFrictionAngle;
- maxFs = std::max((Real) 0, maxFs);
- if (Fs > maxFs) {//Plasticity condition on shear force
- if (!phys->cohesionBroken) {
- phys->cohesionBroken=1;
- phys->shearAdhesion=0;
- maxFs = max((Real) 0, Fn*phys->tangensOfFrictionAngle);
- }
- maxFs = maxFs / Fs;
- shearForce *= maxFs;
- }
-
- //rotational moment are only applied if the cohesion is not broken.
- /// Twist /// the twist law is driven by twist displacement ("getTwist()").
- if(!phys->cohesionBroken){
- Real twist = geom->getTwist() - phys->twp;
- Real twistM=twist*phys->ktw; //elastic twist moment.
- bool sgnChanged=0; //whether the twist moment just passed the equilibrium state.
- if(!contact->isFresh(scene) && phys->moment_twist.dot(twistM*geom->normal)<0)sgnChanged=1;
- if(std::abs(twist)>phys->maxTwist){
- phys->cohesionBroken=1;
- twistM=0;
- }
- else{
- if(std::abs(twistM)>phys->maxElTw || phys->onPlastTw){ //plastic deformation.
- phys->onPlastTw=1;
- if(std::abs(phys->maxCrpRchdTw[0])>std::abs(twist)){ //unloading/reloading
- twistM = phys->kTwUnld*(twist-phys->maxCrpRchdTw[0])+phys->maxCrpRchdTw[1];
- }
- else{//creep loading.
- int sign = twist<0?-1:1;
- twistM = sign*phys->maxElTw+phys->kTwCrp*(twist-sign*phys->maxElTw/phys->ktw); //creep
- phys->maxCrpRchdTw[0]=twist; //new maximum reached
- phys->maxCrpRchdTw[1]=twistM;
- }
- if(sgnChanged){
- phys->maxElTw=0;
- phys->twp=geom->getTwist();
- phys->maxCrpRchdTw[0]=0;
- }
- }
- else{ //elasticity
- phys->maxCrpRchdTw[0]=twist;
- phys->maxCrpRchdTw[1]=twistM;
- }
- }
- phys->moment_twist = twistM * geom->normal;
- }
- else phys->moment_twist=Vector3r(0,0,0);
-
- /// Bending /// incremental form.
- if(!phys->cohesionBroken){
- Vector3r bendM = phys->moment_bending;
- Vector3r relAngVel = geom->getRelAngVel(de1,de2,dt);
- Vector3r relRotBend = (relAngVel - geom->normal.dot(relAngVel)*geom->normal)*dt; // relative rotation due to rolling behaviour
- bendM = geom->rotate(phys->moment_bending); // rotate moment vector (updated)
- phys->pureCreep=geom->rotate(phys->pureCreep); // pure creep is updated to compute the damage.
- Vector3r bendM_elast = bendM-phys->kr*relRotBend;
- if(bendM_elast.norm()>phys->maxElB || phys->onPlastB){ // plastic behavior
- phys->onPlastB=1;
- bendM=bendM-phys->kDam*relRotBend; //trial bending
- if(bendM.norm()<phys->moment_bending.norm()){ // if bending decreased, we are unloading ...
- bendM = bendM+phys->kDam*relRotBend-phys->kRUnld*relRotBend; // ... so undo bendM and apply unload coefficient.
- Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend; // trial pure creep.
- phys->pureCreep = newPureCreep.norm()<phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while unloading, pure creep must decrease.
- phys->kDam=phys->kr+(phys->kRCrp-phys->kr)*(phys->maxCrpRchdB.norm()-phys->maxElB)/(phys->maxBendMom-phys->maxElB); // compute the damage coefficient.
- }
- else{ // bending increased, so we are loading (bendM has to be unchanged).
- Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend;
- phys->pureCreep = newPureCreep.norm()>phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while loading, pure creep must increase.
- if(phys->pureCreep.norm()<bendM.norm()) bendM=phys->pureCreep; // bending moment can't be greather than pure creep.
- if(phys->pureCreep.norm()>phys->maxCrpRchdB.norm()) phys->maxCrpRchdB=phys->pureCreep; // maxCrpRchdB must follow the maximum of pure creep.
- if(phys->pureCreep.norm()>phys->maxBendMom){
- phys->cohesionBroken=1;
- bendM=bendM_elast=Vector3r(0,0,0);
- }
- }
- phys->moment_bending=bendM;
- }
- else{//elasticity
- phys->pureCreep=phys->moment_bending=phys->maxCrpRchdB=bendM_elast;
- phys->kDam=phys->kRCrp;
- }
- }
- phys->shearForce=shearForce;
- phys->normalForce=-Fn*geom->normal;
- applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
- scene->forces.addTorque(id1,-phys->moment_bending-phys->moment_twist);
- scene->forces.addTorque(id2,phys->moment_bending+phys->moment_twist);
-}
-
-YADE_PLUGIN((Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment));
-
-
-
-
-
-
-
=== removed file 'pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp'
--- pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 2014-05-26 09:09:48 +0000
+++ pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 1970-01-01 00:00:00 +0000
@@ -1,26 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
-* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#pragma once
-#include "CohesiveFrictionalContactLaw.hpp"
-#include "InelastCohFrictPhys.hpp"
-
-class Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment: public LawFunctor{
- public:
- Real normElastEnergy();
- Real shearElastEnergy();
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment,LawFunctor,"This law is currently under developpement. Final version and documentation will come before the end of 2014.",
- ,,
- .def("normElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
- .def("shearElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
- );
- FUNCTOR2D(ScGeom6D,InelastCohFrictPhys);
- DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment);
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.cpp 2014-07-14 20:08:33 +0000
@@ -16,7 +16,7 @@
#include <yade/pkg/dem/ScGeom.hpp>
#include <yade/pkg/dem/CapillaryPhys.hpp>
-#include <yade/pkg/dem/Ip2_FrictMat_FrictMat_MindlinCapillaryPhys.hpp>
+#include <yade/pkg/dem/HertzMindlin.hpp>
#include <yade/core/Omega.hpp>
#include <yade/core/Scene.hpp>
#include <yade/lib/base/Math.hpp>
=== modified file 'pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp'
--- pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/Law2_ScGeom_CapillaryPhys_Capillarity.hpp 2014-07-14 20:08:34 +0000
@@ -11,8 +11,6 @@
#pragma once
#include <yade/core/GlobalEngine.hpp>
-#include <boost/tuple/tuple.hpp>
-#include <utility>
/**
This law allows one to take into account capillary forces/effects between spheres coming from the presence of interparticular liquid bridges (menisci).
=== modified file 'pkg/dem/MicroMacroAnalyser.hpp'
--- pkg/dem/MicroMacroAnalyser.hpp 2014-03-21 18:47:45 +0000
+++ pkg/dem/MicroMacroAnalyser.hpp 2014-07-14 20:08:34 +0000
@@ -10,10 +10,6 @@
#include<yade/core/GlobalEngine.hpp>
#include<yade/lib/triangulation/KinematicLocalisationAnalyser.hpp>
-#include <set>
-#include <boost/tuple/tuple.hpp>
-#include <string>
-#include <fstream>
/*! \brief compute fabric tensor, local porosity, local deformation, and other micromechanicaly defined quantities based on triangulation/tesselation of the packing.
=== removed file 'pkg/dem/NormalInelasticMat.cpp'
--- pkg/dem/NormalInelasticMat.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NormalInelasticMat.cpp 1970-01-01 00:00:00 +0000
@@ -1,12 +0,0 @@
-
-#include "NormalInelasticMat.hpp"
-
-
-NormalInelasticMat::~NormalInelasticMat()
-{
-}
-
-YADE_PLUGIN((NormalInelasticMat));
-
-
-
=== removed file 'pkg/dem/NormalInelasticMat.hpp'
--- pkg/dem/NormalInelasticMat.hpp 2011-05-04 12:37:24 +0000
+++ pkg/dem/NormalInelasticMat.hpp 1970-01-01 00:00:00 +0000
@@ -1,32 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2010 by Jerome Duriez <jerome.duriez@xxxxxxxxxxx> *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-
-#pragma once
-
-
-#include<yade/pkg/common/ElastMat.hpp>
-
-
-class NormalInelasticMat : public FrictMat
-{
- public :
- virtual ~NormalInelasticMat ();
-
-/// Serialization
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticMat,FrictMat,"Material class for particles whose contact obey to a normal inelasticity (governed by this :yref:`coeff_dech<NormalInelasticMat::coeff_dech>`).",
- ((Real,coeff_dech,1.0,,"=kn(unload) / kn(load)"))
- ,
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(NormalInelasticMat,FrictMat);
-};
-
-REGISTER_SERIALIZABLE(NormalInelasticMat);
-
-
=== added file 'pkg/dem/NormalInelasticPM.cpp'
--- pkg/dem/NormalInelasticPM.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/NormalInelasticPM.cpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,198 @@
+/*************************************************************************
+* Copyright (C) 2008 by Jérôme DURIEZ *
+* duriez@xxxxxxxxxxxxxxx *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#include "NormalInelasticPM.hpp"
+#include <yade/pkg/dem/ScGeom.hpp>
+#include <yade/core/Omega.hpp>
+#include <yade/core/Scene.hpp>
+
+YADE_PLUGIN((NormalInelasticMat)(NormalInelasticityPhys)(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity)(Ip2_2xNormalInelasticMat_NormalInelasticityPhys));
+
+
+void Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
+{
+ int id1 = contact->getId1();
+ int id2 = contact->getId2();
+// cout << "contact entre " << id1 << " et " << id2;
+
+ NormalInelasticMat* Mat1 = static_cast<NormalInelasticMat*>(Body::byId(id1,scene)->material.get());
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*>(iG.get());
+ NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (iP.get());
+
+ Vector3r& shearForce = currentContactPhysics->shearForce;
+
+ if (contact->isFresh(scene))
+ {
+ shearForce = Vector3r::Zero();
+ currentContactPhysics->previousun=0.0;
+ currentContactPhysics->previousFn=0.0;
+ currentContactPhysics->unMax=0.0;
+ }
+
+ un = geom->penetrationDepth; // >0 for real penetration
+
+// Check if there is a real overlap or not. The Ig2... seems to let exist interactions with negative un (= no overlap). Such interactions seem then to have to be deleted here.
+ if ( un < 0 )
+ {
+ scene->interactions->requestErase(contact);// this, among other things, resets the interaction : geometry and physics variables (as forces, ...) are reset to defaut values
+ return;
+ }
+
+
+
+// ******** Computation of normal Force : depends of the history ******* //
+
+// cout << " Dans Law2 valeur de kn : " << currentContactPhysics->kn << endl;
+// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << currentContactPhysics->previousun << " et previousFn =" << currentContactPhysics->previousFn << endl;
+
+ if(un >= currentContactPhysics->unMax) // case of virgin load : on the "principal line" (limit state of the (un,Fn) space)
+ {
+ Fn = currentContactPhysics->knLower*un;
+ currentContactPhysics->unMax = std::abs(un);
+// cout << "je suis dans le calcul normal " << endl;
+ }
+ else// a priori then we need a greater stifness. False in the case when we are already on the limit state, but a correction below will then do the job
+ {
+ currentContactPhysics->kn = currentContactPhysics->knLower* Mat1->coeff_dech;
+ // Incremental computation of the new normal force :
+ Fn = currentContactPhysics->previousFn + currentContactPhysics->kn * (un-currentContactPhysics->previousun);
+// cout << "je suis dans l'autre calcul" << endl;
+ if(std::abs(Fn) > std::abs(currentContactPhysics->knLower * un)) // check if the limit state is not violated
+ {
+ Fn = currentContactPhysics->knLower*un;
+// cout << "j'etais dans l'autre calcul mais j'ai corrige pour ne pas depasser la limite" << endl;
+ }
+ if(Fn < 0.0 ) // check to stay >=0
+ {
+ Fn = 0;
+// cout << "j'ai corrige pour ne pas etre negatif" << endl;
+ }
+ }
+ currentContactPhysics->normalForce = Fn*geom->normal;
+// cout << "Fn appliquee " << Fn << endl << endl;
+
+ // actualisation :
+ currentContactPhysics->previousFn = Fn;
+ currentContactPhysics->previousun = un;
+
+// *** End of computation of normal force *** //
+
+
+
+// ******** Tangential force ******* //
+ if ( un < 0 )
+ { // BREAK due to tension
+ scene->interactions->requestErase(contact); return;
+ // probably not useful anymore
+ currentContactPhysics->normalForce = Vector3r::Zero();
+ currentContactPhysics->shearForce = Vector3r::Zero();
+ }
+ else
+ {
+ // Correction of previous shear force to take into account the change in normal:
+ shearForce = geom->rotate(currentContactPhysics->shearForce);
+ // Update of shear force corresponding to shear displacement increment:
+ shearForce -= currentContactPhysics->ks * geom->shearIncrement();
+
+ Fs = currentContactPhysics->shearForce.norm();
+ maxFs = std::max((Real) 0,Fn*currentContactPhysics->tangensOfFrictionAngle);
+
+ if ( Fs > maxFs )
+ {
+ maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
+
+ maxFs = maxFs / Fs;
+ if (maxFs>1)
+ cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!" << endl;
+ shearForce *= maxFs;
+ if (Fn<0) currentContactPhysics->normalForce = Vector3r::Zero();
+ }
+
+
+ f = currentContactPhysics->normalForce + shearForce;
+ scene->forces.addForce (id1,-f);
+ scene->forces.addForce (id2, f);
+ scene->forces.addTorque(id1,-(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(f));
+ scene->forces.addTorque(id2, -(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(f));
+// *** End of computation of tangential force *** //
+
+
+// ******** Moment law ******* //
+
+ if(momentRotationLaw)
+ {
+ currentContactPhysics->moment_twist = (geom->getTwist()*currentContactPhysics->kr)*geom->normal ;
+ currentContactPhysics->moment_bending = geom->getBending() * currentContactPhysics->kr;
+
+ moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
+
+// Limitation by plastic threshold of this part of the moment caused by relative twist and bending
+ if (!momentAlwaysElastic)
+ {
+ Real normeMomentMax = currentContactPhysics->forMaxMoment * std::abs(Fn);
+ if(moment.norm()>normeMomentMax)
+ {
+ moment *= normeMomentMax/moment.norm();
+ }
+ }
+ scene->forces.addTorque(id1,-moment);
+ scene->forces.addTorque(id2, moment);
+ }
+// ******** Moment law END ******* //
+ }
+}
+
+
+void Ip2_2xNormalInelasticMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // NormalInelasticMat
+ , const shared_ptr<Material>& b2 // NormalInelasticMat
+ , const shared_ptr<Interaction>& interaction)
+{
+ NormalInelasticMat* sdec1 = static_cast<NormalInelasticMat*>(b1.get());
+ NormalInelasticMat* sdec2 = static_cast<NormalInelasticMat*>(b2.get());
+ ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
+
+
+ if(geom) // so it is ScGeom - NON PERMANENT LINK
+ {
+ if(!interaction->phys)
+ {
+//std::cerr << " isNew, id1: " << interaction->getId1() << " id2: " << interaction->getId2() << "\n";
+ interaction->phys = shared_ptr<NormalInelasticityPhys>(new NormalInelasticityPhys());
+ NormalInelasticityPhys* contactPhysics = YADE_CAST<NormalInelasticityPhys*>(interaction->phys.get());
+
+ Real Ea = sdec1->young;
+ Real Eb = sdec2->young;
+ Real Va = sdec1->poisson;
+ Real Vb = sdec2->poisson;
+ Real Ra = geom->radius1;
+ Real Rb = geom->radius2;
+ Real fa = sdec1->frictionAngle;
+ Real fb = sdec2->frictionAngle;
+
+ Real Kn = 2.0*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);//harmonic average of two stiffnesses
+
+ Real Ks = 2.0*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
+
+ // Jean-Patrick Plassiard, Noura Belheine, Frederic Victor Donze, "A Spherical Discrete Element Model: calibration procedure and incremental response", DOI: 10.1007/s10035-009-0130-x
+
+ Real Kr = betaR*std::pow((Ra+Rb)/2.0,2)*Ks;
+
+ contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
+ contactPhysics->forMaxMoment = 1.0*(Ra+Rb)/2.0; // 1.0 corresponding to ethaR which I don't know exactly where to define as a parameter...
+
+ // Lot of suppress here around (>) r2276.
+ contactPhysics->knLower = Kn;
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
+ contactPhysics->kr = Kr;
+ }
+
+ }
+
+};
+
=== added file 'pkg/dem/NormalInelasticPM.hpp'
--- pkg/dem/NormalInelasticPM.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/NormalInelasticPM.hpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,114 @@
+/*************************************************************************
+* Copyright (C) 2010 by Jerome Duriez <jerome.duriez@xxxxxxxxxxx> *
+* *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+#pragma once
+
+#include <yade/pkg/common/ElastMat.hpp>
+#include <yade/pkg/dem/FrictPhys.hpp>
+#include <yade/pkg/common/Dispatching.hpp>
+#include <yade/pkg/dem/ScGeom.hpp>
+
+class NormalInelasticMat : public FrictMat
+{
+ public :
+ virtual ~NormalInelasticMat () {};
+
+/// Serialization
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticMat,FrictMat,"Material class for particles whose contact obey to a normal inelasticity (governed by this :yref:`coeff_dech<NormalInelasticMat::coeff_dech>`).",
+ ((Real,coeff_dech,1.0,,"=kn(unload) / kn(load)"))
+ ,
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(NormalInelasticMat,FrictMat);
+};
+
+REGISTER_SERIALIZABLE(NormalInelasticMat);
+
+
+class NormalInelasticityPhys : public FrictPhys
+{
+ public :
+ virtual ~NormalInelasticityPhys() {};
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
+ "Physics (of interaction) for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity` : with inelastic unloadings",
+ ((Real,unMax,0.0,,"the maximum value of penetration depth of the history of this interaction"))
+ ((Real,previousun,0.0,,"the value of this un at the last time step"))
+ ((Real,previousFn,0.0,,"the value of the normal force at the last time step"))
+ ((Real,forMaxMoment,1.0,,"parameter stored for each interaction, and allowing to compute the maximum value of the exchanged torque : TorqueMax= forMaxMoment * NormalForce"))
+ ((Real,kr,0.0,,"the rolling stiffness of the interaction"))
+ ((Real,knLower,0.0,,"the stifness corresponding to a virgin load for example"))
+ // internal attributes
+ ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
+ ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
+ ,
+ createIndex();
+ );
+ REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhys);
+};
+
+REGISTER_SERIALIZABLE(NormalInelasticityPhys);
+
+
+class Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
+{
+ private :
+ Vector3r moment // the part of the contact torque of the interaction due to relative rotations (a first part is due to contact forces)
+ ,f// contact force
+ ;
+ Real Fn // value of normal force in the interaction
+ ,Fs // shear force
+ ,maxFs; // maximum value of shear force according to Coulomb-like criterion
+ Real un; // value of interpenetration in the interaction
+ public :
+ virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+
+ FUNCTOR2D(ScGeom,NormalInelasticityPhys);
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity,
+ LawFunctor,
+ "Contact law used to simulate granular filler in rock joints [Duriez2009a]_, [Duriez2011]_. It includes possibility of cohesion, moment transfer and inelastic compression behaviour (to reproduce the normal inelasticity observed for rock joints, for the latter).\n\n The moment transfer relation corresponds to the adaptation of the work of Plassiard & Belheine (see in [DeghmReport2006]_ for example), which was realized by J. Kozicki, and is now coded in :yref:`ScGeom6D`.\n\n As others :yref:`LawFunctor`, it uses pre-computed data of the interactions (rigidities, friction angles -with their tan()-, orientations of the interactions); this work is done here in :yref:`Ip2_2xNormalInelasticMat_NormalInelasticityPhys`.\n\n To use this you should also use :yref:`NormalInelasticMat` as material type of the bodies.\n\n The effects of this law are illustrated in examples/normalInelasticity-test.py",
+ ((bool,momentRotationLaw,true,,"boolean, true=> computation of a torque (against relative rotation) exchanged between particles"))
+ ((bool,momentAlwaysElastic,false,,"boolean, true=> the part of the contact torque (caused by relative rotations, which is computed only if momentRotationLaw..) is not limited by a plastic threshold"))
+ ,
+ moment=Vector3r::Zero();
+ f=Vector3r::Zero();
+ Fn=0.0;
+ Fs=0.0;
+ maxFs=0.0;
+ un=0.0;
+ );
+
+};
+
+REGISTER_SERIALIZABLE(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity);
+
+/*! \brief The RelationShips for using Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity
+
+In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed.
+WARNING : as in the others Relationships most of the attributes are computed only once : when the interaction is "new"
+ */
+
+class Ip2_2xNormalInelasticMat_NormalInelasticityPhys : public IPhysFunctor
+{
+ public :
+
+ virtual void go( const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+
+ FUNCTOR2D(NormalInelasticMat,NormalInelasticMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_2xNormalInelasticMat_NormalInelasticityPhys,
+ IPhysFunctor,
+ "The RelationShips for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`\n\n In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
+ ((Real,betaR,0.12,,"Parameter for computing the torque-stifness : T-stifness = betaR * Rmoy^2"))
+ );
+};
+
+REGISTER_SERIALIZABLE(Ip2_2xNormalInelasticMat_NormalInelasticityPhys);
+
=== removed file 'pkg/dem/NormalInelasticityLaw.cpp'
--- pkg/dem/NormalInelasticityLaw.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/NormalInelasticityLaw.cpp 1970-01-01 00:00:00 +0000
@@ -1,153 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include<yade/pkg/dem/NormalInelasticityLaw.hpp>
-
-#include<yade/pkg/dem/NormalInelasticMat.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-YADE_PLUGIN((Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity));
-
-
-
-void Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
-{
- int id1 = contact->getId1();
- int id2 = contact->getId2();
-// cout << "contact entre " << id1 << " et " << id2;
-
- NormalInelasticMat* Mat1 = static_cast<NormalInelasticMat*>(Body::byId(id1,scene)->material.get());
- ScGeom6D* geom = YADE_CAST<ScGeom6D*>(iG.get());
- NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (iP.get());
-
- Vector3r& shearForce = currentContactPhysics->shearForce;
-
- if (contact->isFresh(scene))
- {
- shearForce = Vector3r::Zero();
- currentContactPhysics->previousun=0.0;
- currentContactPhysics->previousFn=0.0;
- currentContactPhysics->unMax=0.0;
- }
-
- un = geom->penetrationDepth; // >0 for real penetration
-
-// Check if there is a real overlap or not. The Ig2... seems to let exist interactions with negative un (= no overlap). Such interactions seem then to have to be deleted here.
- if ( un < 0 )
- {
- scene->interactions->requestErase(contact);// this, among other things, resets the interaction : geometry and physics variables (as forces, ...) are reset to defaut values
- return;
- }
-
-
-
-// ******** Computation of normal Force : depends of the history ******* //
-
-// cout << " Dans Law2 valeur de kn : " << currentContactPhysics->kn << endl;
-// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << currentContactPhysics->previousun << " et previousFn =" << currentContactPhysics->previousFn << endl;
-
- if(un >= currentContactPhysics->unMax) // case of virgin load : on the "principal line" (limit state of the (un,Fn) space)
- {
- Fn = currentContactPhysics->knLower*un;
- currentContactPhysics->unMax = std::abs(un);
-// cout << "je suis dans le calcul normal " << endl;
- }
- else// a priori then we need a greater stifness. False in the case when we are already on the limit state, but a correction below will then do the job
- {
- currentContactPhysics->kn = currentContactPhysics->knLower* Mat1->coeff_dech;
- // Incremental computation of the new normal force :
- Fn = currentContactPhysics->previousFn + currentContactPhysics->kn * (un-currentContactPhysics->previousun);
-// cout << "je suis dans l'autre calcul" << endl;
- if(std::abs(Fn) > std::abs(currentContactPhysics->knLower * un)) // check if the limit state is not violated
- {
- Fn = currentContactPhysics->knLower*un;
-// cout << "j'etais dans l'autre calcul mais j'ai corrige pour ne pas depasser la limite" << endl;
- }
- if(Fn < 0.0 ) // check to stay >=0
- {
- Fn = 0;
-// cout << "j'ai corrige pour ne pas etre negatif" << endl;
- }
- }
- currentContactPhysics->normalForce = Fn*geom->normal;
-// cout << "Fn appliquee " << Fn << endl << endl;
-
- // actualisation :
- currentContactPhysics->previousFn = Fn;
- currentContactPhysics->previousun = un;
-
-// *** End of computation of normal force *** //
-
-
-
-// ******** Tangential force ******* //
- if ( un < 0 )
- { // BREAK due to tension
- scene->interactions->requestErase(contact); return;
- // probably not useful anymore
- currentContactPhysics->normalForce = Vector3r::Zero();
- currentContactPhysics->shearForce = Vector3r::Zero();
- }
- else
- {
- // Correction of previous shear force to take into account the change in normal:
- shearForce = geom->rotate(currentContactPhysics->shearForce);
- // Update of shear force corresponding to shear displacement increment:
- shearForce -= currentContactPhysics->ks * geom->shearIncrement();
-
- Fs = currentContactPhysics->shearForce.norm();
- maxFs = std::max((Real) 0,Fn*currentContactPhysics->tangensOfFrictionAngle);
-
- if ( Fs > maxFs )
- {
- maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
-
- maxFs = maxFs / Fs;
- if (maxFs>1)
- cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!" << endl;
- shearForce *= maxFs;
- if (Fn<0) currentContactPhysics->normalForce = Vector3r::Zero();
- }
-
-
- f = currentContactPhysics->normalForce + shearForce;
- scene->forces.addForce (id1,-f);
- scene->forces.addForce (id2, f);
- scene->forces.addTorque(id1,-(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(f));
- scene->forces.addTorque(id2, -(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(f));
-// *** End of computation of tangential force *** //
-
-
-// ******** Moment law ******* //
-
- if(momentRotationLaw)
- {
- currentContactPhysics->moment_twist = (geom->getTwist()*currentContactPhysics->kr)*geom->normal ;
- currentContactPhysics->moment_bending = geom->getBending() * currentContactPhysics->kr;
-
- moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
-
-// Limitation by plastic threshold of this part of the moment caused by relative twist and bending
- if (!momentAlwaysElastic)
- {
- Real normeMomentMax = currentContactPhysics->forMaxMoment * std::abs(Fn);
- if(moment.norm()>normeMomentMax)
- {
- moment *= normeMomentMax/moment.norm();
- }
- }
- scene->forces.addTorque(id1,-moment);
- scene->forces.addTorque(id2, moment);
- }
-// ******** Moment law END ******* //
- }
-}
-
-
-
=== removed file 'pkg/dem/NormalInelasticityLaw.hpp'
--- pkg/dem/NormalInelasticityLaw.hpp 2014-05-28 08:17:25 +0000
+++ pkg/dem/NormalInelasticityLaw.hpp 1970-01-01 00:00:00 +0000
@@ -1,51 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/common/Dispatching.hpp>
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/NormalInelasticityPhys.hpp>
-#include <set>
-#include <boost/tuple/tuple.hpp>
-
-
-
-class Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
-{
- private :
- Vector3r moment // the part of the contact torque of the interaction due to relative rotations (a first part is due to contact forces)
- ,f// contact force
- ;
- Real Fn // value of normal force in the interaction
- ,Fs // shear force
- ,maxFs; // maximum value of shear force according to Coulomb-like criterion
- Real un; // value of interpenetration in the interaction
- public :
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
-
- FUNCTOR2D(ScGeom,NormalInelasticityPhys);
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity,
- LawFunctor,
- "Contact law used to simulate granular filler in rock joints [Duriez2009a]_, [Duriez2011]_. It includes possibility of cohesion, moment transfer and inelastic compression behaviour (to reproduce the normal inelasticity observed for rock joints, for the latter).\n\n The moment transfer relation corresponds to the adaptation of the work of Plassiard & Belheine (see in [DeghmReport2006]_ for example), which was realized by J. Kozicki, and is now coded in :yref:`ScGeom6D`.\n\n As others :yref:`LawFunctor`, it uses pre-computed data of the interactions (rigidities, friction angles -with their tan()-, orientations of the interactions); this work is done here in :yref:`Ip2_2xNormalInelasticMat_NormalInelasticityPhys`.\n\n To use this you should also use :yref:`NormalInelasticMat` as material type of the bodies.\n\n The effects of this law are illustrated in examples/normalInelasticity-test.py",
- ((bool,momentRotationLaw,true,,"boolean, true=> computation of a torque (against relative rotation) exchanged between particles"))
- ((bool,momentAlwaysElastic,false,,"boolean, true=> the part of the contact torque (caused by relative rotations, which is computed only if momentRotationLaw..) is not limited by a plastic threshold"))
- ,
- moment=Vector3r::Zero();
- f=Vector3r::Zero();
- Fn=0.0;
- Fs=0.0;
- maxFs=0.0;
- un=0.0;
- );
-
-};
-
-REGISTER_SERIALIZABLE(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity);
-
=== removed file 'pkg/dem/NormalInelasticityPhys.cpp'
--- pkg/dem/NormalInelasticityPhys.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NormalInelasticityPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,18 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include "NormalInelasticityPhys.hpp"
-
-
-NormalInelasticityPhys::~NormalInelasticityPhys()
-{
-}
-
-
-YADE_PLUGIN((NormalInelasticityPhys));
-
=== removed file 'pkg/dem/NormalInelasticityPhys.hpp'
--- pkg/dem/NormalInelasticityPhys.hpp 2011-01-24 14:45:35 +0000
+++ pkg/dem/NormalInelasticityPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,38 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* This program is free software; it is licensed under the terms of the *
-* GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include<yade/pkg/dem/FrictPhys.hpp>
-
-
-class NormalInelasticityPhys : public FrictPhys
-{
- public :
- virtual ~NormalInelasticityPhys();
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
- "Physics (of interaction) for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity` : with inelastic unloadings",
- ((Real,unMax,0.0,,"the maximum value of penetration depth of the history of this interaction"))
- ((Real,previousun,0.0,,"the value of this un at the last time step"))
- ((Real,previousFn,0.0,,"the value of the normal force at the last time step"))
- ((Real,forMaxMoment,1.0,,"parameter stored for each interaction, and allowing to compute the maximum value of the exchanged torque : TorqueMax= forMaxMoment * NormalForce"))
- ((Real,kr,0.0,,"the rolling stiffness of the interaction"))
- ((Real,knLower,0.0,,"the stifness corresponding to a virgin load for example"))
- // internal attributes
- ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
- ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
- ,
- createIndex();
- );
- REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhys);
-};
-
-REGISTER_SERIALIZABLE(NormalInelasticityPhys);
-
-
=== modified file 'pkg/dem/Shop_01.cpp'
--- pkg/dem/Shop_01.cpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/Shop_01.cpp 2014-07-14 20:08:34 +0000
@@ -1,9 +1,6 @@
// 2007 © Václav Šmilauer <eudoxos@xxxxxxxx>
#include"Shop.hpp"
-
-#include<boost/filesystem/convenience.hpp>
#include<boost/tokenizer.hpp>
-#include<boost/tuple/tuple.hpp>
#include"yade/core/Scene.hpp"
#include"yade/core/Body.hpp"
@@ -24,7 +21,7 @@
#include"yade/pkg/dem/NewtonIntegrator.hpp"
#include"yade/pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp"
#include"yade/pkg/dem/Ig2_Box_Sphere_ScGeom.hpp"
-#include"yade/pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp"
+#include"yade/pkg/dem/FrictPhys.hpp"
#include"yade/pkg/common/ForceResetter.hpp"
=== modified file 'pkg/dem/Shop_02.cpp'
--- pkg/dem/Shop_02.cpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/Shop_02.cpp 2014-07-14 20:08:34 +0000
@@ -1,12 +1,6 @@
// 2007 © Václav Šmilauer <eudoxos@xxxxxxxx>
#include"Shop.hpp"
-#include<limits>
-
-#include<boost/filesystem/convenience.hpp>
-#include<boost/tokenizer.hpp>
-#include<boost/tuple/tuple.hpp>
-
#include<yade/core/Scene.hpp>
#include<yade/core/Body.hpp>
#include<yade/core/Interaction.hpp>
@@ -26,7 +20,7 @@
#include<yade/pkg/dem/NewtonIntegrator.hpp>
#include<yade/pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp>
#include<yade/pkg/dem/Ig2_Box_Sphere_ScGeom.hpp>
-#include<yade/pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp>
+#include<yade/pkg/dem/FrictPhys.hpp>
#include<yade/pkg/common/ForceResetter.hpp>
=== modified file 'pkg/dem/SimpleShear.cpp'
--- pkg/dem/SimpleShear.cpp 2014-07-02 16:18:24 +0000
+++ pkg/dem/SimpleShear.cpp 2014-07-14 20:08:34 +0000
@@ -11,9 +11,7 @@
#include "SimpleShear.hpp"
-#include <yade/pkg/dem/NormalInelasticMat.hpp>
-#include<yade/pkg/dem/NormalInelasticityLaw.hpp>
-#include <yade/pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp>
+#include <yade/pkg/dem/NormalInelasticPM.hpp>
#include<yade/pkg/dem/GlobalStiffnessTimeStepper.hpp>
#include<yade/pkg/common/Aabb.hpp>
@@ -27,7 +25,7 @@
#include<yade/pkg/dem/NewtonIntegrator.hpp>
#include<yade/pkg/common/GravityEngines.hpp>
-#include<yade/pkg/dem/KinemCTDEngine.hpp>
+#include<yade/pkg/dem/KinemC__Engine.hpp>
#include<yade/pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp>
#include<yade/pkg/dem/Ig2_Box_Sphere_ScGeom.hpp>
@@ -38,7 +36,6 @@
#include<yade/pkg/common/Box.hpp>
#include<yade/pkg/common/Sphere.hpp>
-#include <boost/filesystem/convenience.hpp>
#include <utility>
YADE_PLUGIN((SimpleShear))
=== modified file 'pkg/dem/TriaxialTest.cpp'
--- pkg/dem/TriaxialTest.cpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/TriaxialTest.cpp 2014-07-14 20:08:34 +0000
@@ -8,7 +8,7 @@
#include<yade/pkg/dem/ElasticContactLaw.hpp>
-#include<yade/pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp>
+#include<yade/pkg/dem/FrictPhys.hpp>
#include<yade/pkg/dem/GlobalStiffnessTimeStepper.hpp>
#include<yade/pkg/common/ElastMat.hpp>
#include<yade/pkg/dem/TriaxialStressController.hpp>
@@ -34,13 +34,11 @@
#include<yade/pkg/dem/Ig2_Sphere_Sphere_ScGeom.hpp>
#include<yade/pkg/dem/Ig2_Box_Sphere_ScGeom.hpp>
#include<yade/pkg/dem/Ig2_Facet_Sphere_ScGeom.hpp>
-#include<yade/pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.hpp>
#include<yade/pkg/common/Bo1_Sphere_Aabb.hpp>
#include<yade/pkg/common/Bo1_Box_Aabb.hpp>
#include<yade/pkg/common/Bo1_Facet_Aabb.hpp>
#include<yade/pkg/common/Wall.hpp>
-#include <boost/filesystem/convenience.hpp>
#include <boost/numeric/conversion/bounds.hpp>
#include <boost/limits.hpp>
// random
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/VTKRecorder.cpp 2014-07-16 15:50:43 +0000
@@ -7,7 +7,7 @@
#include<vtkPointData.h>
#include<vtkCellData.h>
#include<vtkSmartPointer.h>
-#include<vtkFloatArray.h>
+#include<vtkDoubleArray.h>
#include<vtkUnstructuredGrid.h>
#include<vtkPolyData.h>
#include<vtkXMLUnstructuredGridWriter.h>
@@ -79,7 +79,8 @@
else if(rec=="cracks") recActive[REC_CRACKS]=true;
else if(rec=="pericell" && scene->isPeriodic) recActive[REC_PERICELL]=true;
else if(rec=="liquidcontrol") recActive[REC_LIQ]=true;
- else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, jcfpm, cracks, pericell). Ignored.");
+ else if(rec=="bstresses") recActive[REC_BSTRESS]=true;
+ else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, jcfpm, cracks, pericell, liquidcontrol, bstresses). Ignored.");
}
// cpm needs interactions
if(recActive[REC_CPM]) recActive[REC_INTR]=true;
@@ -98,147 +99,171 @@
vtkSmartPointer<vtkPoints> spheresPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> spheresCells = vtkSmartPointer<vtkCellArray>::New();
- vtkSmartPointer<vtkFloatArray> radii = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> radii = vtkSmartPointer<vtkDoubleArray>::New();
radii->SetNumberOfComponents(1);
radii->SetName("radii");
- vtkSmartPointer<vtkFloatArray> spheresMass = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresSigI = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresSigI->SetNumberOfComponents(1);
+ spheresSigI->SetName("sigI");
+
+ vtkSmartPointer<vtkDoubleArray> spheresSigII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresSigII->SetNumberOfComponents(1);
+ spheresSigII->SetName("sigII");
+
+ vtkSmartPointer<vtkDoubleArray> spheresSigIII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresSigIII->SetNumberOfComponents(1);
+ spheresSigIII->SetName("sigIII");
+
+ vtkSmartPointer<vtkDoubleArray> spheresDirI = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresDirI->SetNumberOfComponents(3);
+ spheresDirI->SetName("dirI");
+
+ vtkSmartPointer<vtkDoubleArray> spheresDirII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresDirII->SetNumberOfComponents(3);
+ spheresDirII->SetName("dirII");
+
+ vtkSmartPointer<vtkDoubleArray> spheresDirIII = vtkSmartPointer<vtkDoubleArray>::New();
+ spheresDirIII->SetNumberOfComponents(3);
+ spheresDirIII->SetName("dirIII");
+
+ vtkSmartPointer<vtkDoubleArray> spheresMass = vtkSmartPointer<vtkDoubleArray>::New();
spheresMass->SetNumberOfComponents(1);
spheresMass->SetName("mass");
- vtkSmartPointer<vtkFloatArray> spheresId = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresId = vtkSmartPointer<vtkDoubleArray>::New();
spheresId->SetNumberOfComponents(1);
spheresId->SetName("id");
#ifdef YADE_SPH
- vtkSmartPointer<vtkFloatArray> spheresCsSPH = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresCsSPH = vtkSmartPointer<vtkDoubleArray>::New();
spheresCsSPH->SetNumberOfComponents(1);
spheresCsSPH->SetName("SPH_Cs");
- vtkSmartPointer<vtkFloatArray> spheresRhoSPH = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresRhoSPH = vtkSmartPointer<vtkDoubleArray>::New();
spheresRhoSPH->SetNumberOfComponents(1);
spheresRhoSPH->SetName("SPH_Rho");
- vtkSmartPointer<vtkFloatArray> spheresPressSPH = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresPressSPH = vtkSmartPointer<vtkDoubleArray>::New();
spheresPressSPH->SetNumberOfComponents(1);
spheresPressSPH->SetName("SPH_Press");
- vtkSmartPointer<vtkFloatArray> spheresCoordNumbSPH = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresCoordNumbSPH = vtkSmartPointer<vtkDoubleArray>::New();
spheresCoordNumbSPH->SetNumberOfComponents(1);
spheresCoordNumbSPH->SetName("SPH_Neigh");
#endif
#ifdef YADE_LIQMIGRATION
- vtkSmartPointer<vtkFloatArray> spheresLiqVol = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresLiqVol = vtkSmartPointer<vtkDoubleArray>::New();
spheresLiqVol->SetNumberOfComponents(1);
spheresLiqVol->SetName("Liq_Vol");
- vtkSmartPointer<vtkFloatArray> spheresLiqVolIter = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresLiqVolIter = vtkSmartPointer<vtkDoubleArray>::New();
spheresLiqVolIter->SetNumberOfComponents(1);
spheresLiqVolIter->SetName("Liq_VolIter");
- vtkSmartPointer<vtkFloatArray> spheresLiqVolTotal = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresLiqVolTotal = vtkSmartPointer<vtkDoubleArray>::New();
spheresLiqVolTotal->SetNumberOfComponents(1);
spheresLiqVolTotal->SetName("Liq_VolTotal");
#endif
- vtkSmartPointer<vtkFloatArray> spheresMask = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresMask = vtkSmartPointer<vtkDoubleArray>::New();
spheresMask->SetNumberOfComponents(1);
spheresMask->SetName("mask");
- vtkSmartPointer<vtkFloatArray> clumpId = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> clumpId = vtkSmartPointer<vtkDoubleArray>::New();
clumpId->SetNumberOfComponents(1);
clumpId->SetName("clumpId");
- vtkSmartPointer<vtkFloatArray> spheresColors = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresColors = vtkSmartPointer<vtkDoubleArray>::New();
spheresColors->SetNumberOfComponents(3);
spheresColors->SetName("color");
- vtkSmartPointer<vtkFloatArray> spheresLinVelVec = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresLinVelVec = vtkSmartPointer<vtkDoubleArray>::New();
spheresLinVelVec->SetNumberOfComponents(3);
spheresLinVelVec->SetName("linVelVec"); //Linear velocity in Vector3 form
- vtkSmartPointer<vtkFloatArray> spheresLinVelLen = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresLinVelLen = vtkSmartPointer<vtkDoubleArray>::New();
spheresLinVelLen->SetNumberOfComponents(1);
spheresLinVelLen->SetName("linVelLen"); //Length (magnitude) of linear velocity
- vtkSmartPointer<vtkFloatArray> spheresAngVelVec = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresAngVelVec = vtkSmartPointer<vtkDoubleArray>::New();
spheresAngVelVec->SetNumberOfComponents(3);
spheresAngVelVec->SetName("angVelVec"); //Angular velocity in Vector3 form
- vtkSmartPointer<vtkFloatArray> spheresAngVelLen = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresAngVelLen = vtkSmartPointer<vtkDoubleArray>::New();
spheresAngVelLen->SetNumberOfComponents(1);
spheresAngVelLen->SetName("angVelLen"); //Length (magnitude) of angular velocity
- vtkSmartPointer<vtkFloatArray> spheresNormalStressVec = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresNormalStressVec = vtkSmartPointer<vtkDoubleArray>::New();
spheresNormalStressVec->SetNumberOfComponents(3);
spheresNormalStressVec->SetName("normalStress");
- vtkSmartPointer<vtkFloatArray> spheresShearStressVec = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresShearStressVec = vtkSmartPointer<vtkDoubleArray>::New();
spheresShearStressVec->SetNumberOfComponents(3);
spheresShearStressVec->SetName("shearStress");
- vtkSmartPointer<vtkFloatArray> spheresNormalStressNorm = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresNormalStressNorm = vtkSmartPointer<vtkDoubleArray>::New();
spheresNormalStressNorm->SetNumberOfComponents(1);
spheresNormalStressNorm->SetName("normalStressNorm");
- vtkSmartPointer<vtkFloatArray> spheresMaterialId = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> spheresMaterialId = vtkSmartPointer<vtkDoubleArray>::New();
spheresMaterialId->SetNumberOfComponents(1);
spheresMaterialId->SetName("materialId");
// facets
vtkSmartPointer<vtkPoints> facetsPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> facetsCells = vtkSmartPointer<vtkCellArray>::New();
- vtkSmartPointer<vtkFloatArray> facetsColors = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> facetsColors = vtkSmartPointer<vtkDoubleArray>::New();
facetsColors->SetNumberOfComponents(3);
facetsColors->SetName("color");
- vtkSmartPointer<vtkFloatArray> facetsForceVec = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> facetsForceVec = vtkSmartPointer<vtkDoubleArray>::New();
facetsForceVec->SetNumberOfComponents(3);
facetsForceVec->SetName("stressVec");
- vtkSmartPointer<vtkFloatArray> facetsForceLen = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> facetsForceLen = vtkSmartPointer<vtkDoubleArray>::New();
facetsForceLen->SetNumberOfComponents(1);
facetsForceLen->SetName("stressLen");
- vtkSmartPointer<vtkFloatArray> facetsMaterialId = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> facetsMaterialId = vtkSmartPointer<vtkDoubleArray>::New();
facetsMaterialId->SetNumberOfComponents(1);
facetsMaterialId->SetName("materialId");
- vtkSmartPointer<vtkFloatArray> facetsMask = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> facetsMask = vtkSmartPointer<vtkDoubleArray>::New();
facetsMask->SetNumberOfComponents(1);
facetsMask->SetName("mask");
// boxes
vtkSmartPointer<vtkPoints> boxesPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> boxesCells = vtkSmartPointer<vtkCellArray>::New();
- vtkSmartPointer<vtkFloatArray> boxesColors = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> boxesColors = vtkSmartPointer<vtkDoubleArray>::New();
boxesColors->SetNumberOfComponents(3);
boxesColors->SetName("color");
- vtkSmartPointer<vtkFloatArray> boxesForceVec = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> boxesForceVec = vtkSmartPointer<vtkDoubleArray>::New();
boxesForceVec->SetNumberOfComponents(3);
boxesForceVec->SetName("stressVec");
- vtkSmartPointer<vtkFloatArray> boxesForceLen = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> boxesForceLen = vtkSmartPointer<vtkDoubleArray>::New();
boxesForceLen->SetNumberOfComponents(1);
boxesForceLen->SetName("stressLen");
- vtkSmartPointer<vtkFloatArray> boxesMaterialId = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> boxesMaterialId = vtkSmartPointer<vtkDoubleArray>::New();
boxesMaterialId->SetNumberOfComponents(1);
boxesMaterialId->SetName("materialId");
- vtkSmartPointer<vtkFloatArray> boxesMask = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> boxesMask = vtkSmartPointer<vtkDoubleArray>::New();
boxesMask->SetNumberOfComponents(1);
boxesMask->SetName("mask");
// interactions
vtkSmartPointer<vtkPoints> intrBodyPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> intrCells = vtkSmartPointer<vtkCellArray>::New();
- vtkSmartPointer<vtkFloatArray> intrForceN = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> intrForceN = vtkSmartPointer<vtkDoubleArray>::New();
intrForceN->SetNumberOfComponents(1);
intrForceN->SetName("forceN");
- vtkSmartPointer<vtkFloatArray> intrAbsForceT = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> intrAbsForceT = vtkSmartPointer<vtkDoubleArray>::New();
intrAbsForceT->SetNumberOfComponents(3);
intrAbsForceT->SetName("absForceT");
@@ -248,73 +273,58 @@
// extras for CPM
if(recActive[REC_CPM]){ CpmStateUpdater csu; csu.update(scene); }
- vtkSmartPointer<vtkFloatArray> cpmDamage = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> cpmDamage = vtkSmartPointer<vtkDoubleArray>::New();
cpmDamage->SetNumberOfComponents(1);
cpmDamage->SetName("cpmDamage");
- vtkSmartPointer<vtkFloatArray> cpmStress = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> cpmStress = vtkSmartPointer<vtkDoubleArray>::New();
cpmStress->SetNumberOfComponents(9);
cpmStress->SetName("cpmStress");
// extras for JCFpm
- vtkSmartPointer<vtkFloatArray> damage = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> damage = vtkSmartPointer<vtkDoubleArray>::New();
damage->SetNumberOfComponents(1);;
damage->SetName("damage");
- vtkSmartPointer<vtkFloatArray> damageRel = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> damageRel = vtkSmartPointer<vtkDoubleArray>::New();
damageRel->SetNumberOfComponents(1);;
damageRel->SetName("damageRel");
- vtkSmartPointer<vtkFloatArray> intrIsCohesive = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> intrIsCohesive = vtkSmartPointer<vtkDoubleArray>::New();
intrIsCohesive->SetNumberOfComponents(1);
intrIsCohesive->SetName("isCohesive");
- vtkSmartPointer<vtkFloatArray> intrIsOnJoint = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> intrIsOnJoint = vtkSmartPointer<vtkDoubleArray>::New();
intrIsOnJoint->SetNumberOfComponents(1);
intrIsOnJoint->SetName("isOnJoint");
// extras for cracks
vtkSmartPointer<vtkPoints> crackPos = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> crackCells = vtkSmartPointer<vtkCellArray>::New();
- vtkSmartPointer<vtkFloatArray> iter = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> iter = vtkSmartPointer<vtkDoubleArray>::New();
iter->SetNumberOfComponents(1);
iter->SetName("iter");
- vtkSmartPointer<vtkFloatArray> crackType = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> crackType = vtkSmartPointer<vtkDoubleArray>::New();
crackType->SetNumberOfComponents(1);
crackType->SetName("crackType");
- vtkSmartPointer<vtkFloatArray> crackSize = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> crackSize = vtkSmartPointer<vtkDoubleArray>::New();
crackSize->SetNumberOfComponents(1);
crackSize->SetName("crackSize");
- vtkSmartPointer<vtkFloatArray> crackNorm = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> crackNorm = vtkSmartPointer<vtkDoubleArray>::New();
crackNorm->SetNumberOfComponents(3);
crackNorm->SetName("crackNorm");
#ifdef YADE_LIQMIGRATION
- vtkSmartPointer<vtkFloatArray> liqVol = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> liqVol = vtkSmartPointer<vtkDoubleArray>::New();
liqVol->SetNumberOfComponents(1);
liqVol->SetName("liqVol");
- vtkSmartPointer<vtkFloatArray> liqVolNorm = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> liqVolNorm = vtkSmartPointer<vtkDoubleArray>::New();
liqVolNorm->SetNumberOfComponents(1);
liqVolNorm->SetName("liqVolNorm");
#endif
- // the same for newly created cracks
-// vtkSmartPointer<vtkPoints> crackPosNew = vtkSmartPointer<vtkPoints>::New();
-// vtkSmartPointer<vtkCellArray> crackCellsNew = vtkSmartPointer<vtkCellArray>::New();
-// vtkSmartPointer<vtkFloatArray> iterNew = vtkSmartPointer<vtkFloatArray>::New();
-// iterNew->SetNumberOfComponents(1);
-// iterNew->SetName("iter");
-// vtkSmartPointer<vtkFloatArray> crackTypeNew = vtkSmartPointer<vtkFloatArray>::New();
-// crackTypeNew->SetNumberOfComponents(1);
-// crackTypeNew->SetName("crackType");
-// vtkSmartPointer<vtkFloatArray> crackSizeNew = vtkSmartPointer<vtkFloatArray>::New();
-// crackSizeNew->SetNumberOfComponents(1);
-// crackSizeNew->SetName("crackSize");
-// vtkSmartPointer<vtkFloatArray> crackNormNew = vtkSmartPointer<vtkFloatArray>::New();
-// crackNormNew->SetNumberOfComponents(3);
-// crackNormNew->SetName("crackNorm");
// extras for WireMatPM
- vtkSmartPointer<vtkFloatArray> wpmNormalForce = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> wpmNormalForce = vtkSmartPointer<vtkDoubleArray>::New();
wpmNormalForce->SetNumberOfComponents(1);
wpmNormalForce->SetName("wpmNormalForce");
- vtkSmartPointer<vtkFloatArray> wpmLimitFactor = vtkSmartPointer<vtkFloatArray>::New();
+ vtkSmartPointer<vtkDoubleArray> wpmLimitFactor = vtkSmartPointer<vtkDoubleArray>::New();
wpmLimitFactor->SetNumberOfComponents(1);
wpmLimitFactor->SetName("wpmLimitFactor");
@@ -377,8 +387,8 @@
const NormShearPhys* phys = YADE_CAST<NormShearPhys*>(I->phys.get());
const GenericSpheresContact* geom = YADE_CAST<GenericSpheresContact*>(I->geom.get());
// gives _signed_ scalar of normal force, following the convention used in the respective constitutive law
- float fn=phys->normalForce.dot(geom->normal);
- float fs[3]={ (float) std::abs(phys->shearForce[0]), (float) std::abs(phys->shearForce[1]), (float) std::abs(phys->shearForce[2])};
+ Real fn=phys->normalForce.dot(geom->normal);
+ Real fs[3]={ (Real) std::abs(phys->shearForce[0]), (Real) std::abs(phys->shearForce[1]), (Real) std::abs(phys->shearForce[2])};
// add the value once for each interaction object that we created (might be 2 for the periodic boundary)
for(int i=0; i<numAddValues; i++){
intrAbsForceT->InsertNextTupleValue(fs);
@@ -418,6 +428,13 @@
vector<Shop::bodyState> bodyStates;
if(recActive[REC_STRESS]) Shop::getStressForEachBody(bodyStates);
+
+ vector<Matrix3r> bStresses;
+ if (recActive[REC_BSTRESS])
+ {
+ Shop::getStressLWForEachBody(bStresses);
+ }
+
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
if (!b) continue;
if(mask!=0 && !b->maskCompatible(mask)) continue;
@@ -430,31 +447,76 @@
pid[0] = spheresPos->InsertNextPoint(pos[0], pos[1], pos[2]);
spheresCells->InsertNextCell(1,pid);
radii->InsertNextValue(sphere->radius);
+
+ if (recActive[REC_BSTRESS]) {
+ const Matrix3r& bStress = - bStresses[b->getId()]; // compressive states are negativ for getStressLWForEachBody; I want them as positiv
+ Eigen::SelfAdjointEigenSolver<Matrix3r> solver(bStress); // bStress is probably not symmetric (= self-adjoint for real matrices), but the solver hopefully works (considering only one half of bStress). And, moreover, existence of (real) eigenvalues is not sure for not symmetric bStress..
+ Matrix3r dirAll = solver.eigenvectors();
+ Vector3r eigenVal = solver.eigenvalues(); // cf http://eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointEigenSolver.html#a30caf3c3884a7f4a46b8ec94efd23c5e to be sure that eigenVal[i] * dirAll.col(i) = bStress * dirAll.col(i)
+ int whereSigI(-1), whereSigII(-1), whereSigIII(-1); // all whereSig_i are in [0;2] : whereSigI = 2 => sigI=eigenVal[2]
+ if ( eigenVal[0] > std::max(eigenVal[1],eigenVal[2]) ) {
+ whereSigI = 0;
+ if (eigenVal[1]>eigenVal[2]) {
+ whereSigII=1;
+ whereSigIII=2; }
+ else { //eigenVal[0] > eigenVal[2] >= eigenVal[1]
+ whereSigII = 2;
+ whereSigIII=1; } }
+ else { // max(lambda1,lambda2) >= lambda0 // lambda = eigenVal in the comments
+ if (eigenVal[1]>=eigenVal[2]) {//max(lambda1,lambda2) = lambda1 : lambda 1 >= lambda2
+ whereSigI = 1;
+ if (eigenVal[2]>=eigenVal[0]) {//lambda1 >= lambda2 >= lambda0
+ whereSigII=2;
+ whereSigIII=0; }
+ else { //lambda1 >= lambda0 > lambda2
+ whereSigII=0;
+ whereSigIII=2; } }
+ else { //max(lambda1,lambda2) = lambda2 : lambda2 > lambda1
+ whereSigI = 2;
+ if (eigenVal[1] > eigenVal[0]) {
+ whereSigII = 1;
+ whereSigIII = 0; }
+ else {
+ whereSigIII = 1;
+ whereSigII = 0; } } }
+
+ spheresSigI->InsertNextValue(eigenVal[whereSigI]);
+ spheresSigII->InsertNextValue(eigenVal[whereSigII]);
+ spheresSigIII->InsertNextValue(eigenVal[whereSigIII]);
+ Real dirI[3] { (Real) dirAll(0,whereSigI), (Real) dirAll(1,whereSigI), (Real) dirAll(2,whereSigI) };
+ spheresDirI->InsertNextTupleValue(dirI);
+
+ Real dirII[3] { (Real) dirAll(0,whereSigII), (Real) dirAll(1,whereSigII), (Real) dirAll(2,whereSigII) };
+ spheresDirII->InsertNextTupleValue(dirII);
+
+ Real dirIII[3] { (Real) dirAll(0,whereSigIII), (Real) dirAll(1,whereSigIII), (Real) dirAll(2,whereSigIII) };
+ spheresDirIII->InsertNextTupleValue(dirIII); }
+
if (recActive[REC_ID]) spheresId->InsertNextValue(b->getId());
if (recActive[REC_MASK]) spheresMask->InsertNextValue(GET_MASK(b));
if (recActive[REC_MASS]) spheresMass->InsertNextValue(b->state->mass);
if (recActive[REC_CLUMPID]) clumpId->InsertNextValue(b->clumpId);
if (recActive[REC_COLORS]){
const Vector3r& color = sphere->color;
- float c[3] = { (float) color[0], (float) color[1], (float) color[2]};
+ Real c[3] = { (Real) color[0], (Real) color[1], (Real) color[2]};
spheresColors->InsertNextTupleValue(c);
}
if(recActive[REC_VELOCITY]){
const Vector3r& vel = b->state->vel;
- float v[3] = { (float) vel[0], (float) vel[1], (float) vel[2] };
+ Real v[3] = { (Real) vel[0], (Real) vel[1], (Real) vel[2] };
spheresLinVelVec->InsertNextTupleValue(v);
spheresLinVelLen->InsertNextValue(vel.norm());
const Vector3r& angVel = b->state->angVel;
- float av[3] = { (float) angVel[0], (float) angVel[1], (float) angVel[2] };
+ Real av[3] = { (Real) angVel[0], (Real) angVel[1], (Real) angVel[2] };
spheresAngVelVec->InsertNextTupleValue(av);
spheresAngVelLen->InsertNextValue(angVel.norm());
}
if(recActive[REC_STRESS]){
const Vector3r& stress = bodyStates[b->getId()].normStress;
const Vector3r& shear = bodyStates[b->getId()].shearStress;
- float n[3] = { (float) stress[0], (float) stress[1], (float) stress[2] };
- float s[3] = { (float) shear [0], (float) shear [1], (float) shear [2] };
+ Real n[3] = { (Real) stress[0], (Real) stress[1], (Real) stress[2] };
+ Real s[3] = { (Real) shear [0], (Real) shear [1], (Real) shear [2] };
spheresNormalStressVec->InsertNextTupleValue(n);
spheresShearStressVec->InsertNextTupleValue(s);
spheresNormalStressNorm->InsertNextValue(stress.norm());
@@ -463,8 +525,8 @@
if (recActive[REC_CPM]){
cpmDamage->InsertNextValue(YADE_PTR_CAST<CpmState>(b->state)->normDmg);
const Matrix3r& ss=YADE_PTR_CAST<CpmState>(b->state)->stress;
- //float s[3]={ss[0],ss[1],ss[2]};
- float s[9]={ (float) ss(0,0), (float) ss(0,1), (float) ss(0,2), (float) ss(1,0), (float) ss(1,1), (float) ss(1,2), (float) ss(2,0), (float) ss(2,1), (float) ss(2,2)};
+ //Real s[3]={ss[0],ss[1],ss[2]};
+ Real s[9]={ (Real) ss(0,0), (Real) ss(0,1), (Real) ss(0,2), (Real) ss(1,0), (Real) ss(1,1), (Real) ss(1,2), (Real) ss(2,0), (Real) ss(2,1), (Real) ss(2,2)};
cpmStress->InsertNextTupleValue(s);
}
@@ -506,12 +568,12 @@
facetsCells->InsertNextCell(tri);
if (recActive[REC_COLORS]){
const Vector3r& color = facet->color;
- float c[3] = { (float) color[0], (float) color[1], (float) color[2]};
+ Real c[3] = { (Real) color[0], (Real) color[1], (Real) color[2]};
facetsColors->InsertNextTupleValue(c);
}
if(recActive[REC_STRESS]){
const Vector3r& stress = bodyStates[b->getId()].normStress+bodyStates[b->getId()].shearStress;
- float s[3] = { (float) stress[0], (float) stress[1], (float) stress[2] };
+ Real s[3] = { (Real) stress[0], (Real) stress[1], (Real) stress[2] };
facetsForceVec->InsertNextTupleValue(s);
facetsForceLen->InsertNextValue(stress.norm());
}
@@ -558,12 +620,12 @@
for(int i=0; i<6; i++){
if (recActive[REC_COLORS]){
const Vector3r& color = box->color;
- float c[3] = { (float) color[0], (float) color[1], (float) color[2]};
+ Real c[3] = { (Real) color[0], (Real) color[1], (Real) color[2]};
boxesColors->InsertNextTupleValue(c);
}
if(recActive[REC_STRESS]){
const Vector3r& stress = bodyStates[b->getId()].normStress+bodyStates[b->getId()].shearStress;
- float s[3] = { (float) stress[0], (float) stress[1], (float) stress[2] };
+ Real s[3] = { (Real) stress[0], (Real) stress[1], (Real) stress[2] };
boxesForceVec->InsertNextTupleValue(s);
boxesForceLen->InsertNextValue(stress.norm());
}
@@ -647,7 +709,15 @@
if (recActive[REC_JCFPM]) {
spheresUg->GetPointData()->AddArray(damage);
}
-
+ if (recActive[REC_BSTRESS])
+ {
+ spheresUg->GetPointData()->AddArray(spheresSigI);
+ spheresUg->GetPointData()->AddArray(spheresSigII);
+ spheresUg->GetPointData()->AddArray(spheresSigIII);
+ spheresUg->GetPointData()->AddArray(spheresDirI);
+ spheresUg->GetPointData()->AddArray(spheresDirII);
+ spheresUg->GetPointData()->AddArray(spheresDirIII);
+ }
if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
#ifdef YADE_VTK_MULTIBLOCK
@@ -786,14 +856,12 @@
string fileCracks = "cracks_"+Key+".txt";
std::ifstream file (fileCracks.c_str(),std::ios::in);
vtkSmartPointer<vtkUnstructuredGrid> crackUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
- vtkSmartPointer<vtkUnstructuredGrid> crackUgNew = vtkSmartPointer<vtkUnstructuredGrid>::New();
if(file){
while ( !file.eof() ){
std::string line;
Real i,p0,p1,p2,t,s,n0,n1,n2;
- while ( std::getline(file, line)/* writes into string "line", a line of file "file". To go along diff. lines*/ )
- {
+ while ( std::getline(file, line) {/* writes into string "line", a line of file "file". To go along diff. lines*/ )
file >> i >> p0 >> p1 >> p2 >> t >> s >> n0 >> n1 >> n2;
vtkIdType pid[1];
pid[0] = crackPos->InsertNextPoint(p0, p1, p2);
@@ -801,23 +869,10 @@
crackType->InsertNextValue(t);
crackSize->InsertNextValue(s);
iter->InsertNextValue(i);
- float n[3] = { n0, n1, n2 };
- crackNorm->InsertNextTupleValue(n);
- // Then, taking care only of newly created cracks :
-// if (i > scene->iter - iterPeriod)
-// {
-// pid[0] = crackPosNew->InsertNextPoint(p0, p1, p2);
-// crackCellsNew->InsertNextCell(1,pid);
-// crackTypeNew->InsertNextValue(t);
-// crackSizeNew->InsertNextValue(s);
-// iterNew->InsertNextValue(i);
-// crackNormNew->InsertNextTupleValue(n);
-// }
- }
- }
- file.close();
- }
-
+ Real n[3] = { n0, n1, n2 };
+ crackNorm->InsertNextTupleValue(n); } }
+ file.close(); }
+//
crackUg->SetPoints(crackPos);
crackUg->SetCells(VTK_VERTEX, crackCells);
crackUg->GetPointData()->AddArray(iter);
@@ -835,25 +890,7 @@
#else
writer->SetInput(crackUg);
#endif
- writer->Write();
-
- // Same operations, for newly created cracks :
-// crackUgNew->SetPoints(crackPosNew);
-// crackUgNew->SetCells(VTK_VERTEX, crackCellsNew);
-// crackUgNew->GetPointData()->AddArray(iterNew);
-// crackUgNew->GetPointData()->AddArray(crackTypeNew);
-// crackUgNew->GetPointData()->AddArray(crackSizeNew);
-// crackUgNew->GetPointData()->AddArray(crackNormNew); //same remark about the orientation...
-//
-// fn=fileName+"newcracks."+boost::lexical_cast<string>(scene->iter)+".vtu";
-// writer->SetFileName(fn.c_str());
-// #ifdef YADE_VTK6
-// writer->SetInputData(crackUgNew);
-// #else
-// writer->SetInput(crackUgNew);
-// #endif
-// writer->Write();
- }
+ writer->Write(); }
#ifdef YADE_VTK_MULTIBLOCK
if(multiblock){
=== modified file 'pkg/dem/VTKRecorder.hpp'
--- pkg/dem/VTKRecorder.hpp 2014-06-12 15:11:52 +0000
+++ pkg/dem/VTKRecorder.hpp 2014-07-16 08:34:55 +0000
@@ -11,7 +11,7 @@
class VTKRecorder: public PeriodicEngine {
public:
- enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_JCFPM,REC_CRACKS,REC_WPM,REC_PERICELL,REC_LIQ};
+ enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_JCFPM,REC_CRACKS,REC_WPM,REC_PERICELL,REC_LIQ,REC_BSTRESS};
virtual void action();
void addWallVTK (vtkSmartPointer<vtkQuad>& boxes, vtkSmartPointer<vtkPoints>& boxesPos, Vector3r& W1, Vector3r& W2, Vector3r& W3, Vector3r& W4);
YADE_CLASS_BASE_DOC_ATTRS_CTOR(VTKRecorder,PeriodicEngine,"Engine recording snapshots of simulation into series of \\*.vtu files, readable by VTK-based postprocessing programs such as Paraview. Both bodies (spheres and facets) and interactions can be recorded, with various vector/scalar quantities that are defined on them.\n\n:yref:`PeriodicEngine.initRun` is initialized to ``True`` automatically.",
@@ -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``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n\t``stress``\n\t\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n\t``pericell``\n\t\tSaves the shape of the cell (simulation has to be periodic).\n\n.. admonition:: Specific recorders\n\n\tThe following should only be activated in appropriate cases, otherwise crashes can occur due to violation of type presuppositions.\n\n\t``cpm``\n\t\tSaves data pertaining to the :yref:`concrete model<Law2_ScGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmStress`` (stress on particle); ``intr`` is activated automatically by ``cpm``\n\t``wpm``\n\t\tSaves data pertaining to the :yref:`wire particle model<Law2_ScGeom_WirePhys_WirePM>`: ``wpmForceNFactor`` shows the loading factor for the wire, e.g. normal force divided by threshold normal force.\n\t``jcfpm``\n\t\tSaves data pertaining to the :yref:`rock (smooth)-jointed model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``damage`` is defined by :yref:`JCFpmState.tensBreak` + :yref:`JCFpmState.shearBreak`; ``intr`` is activated automatically by ``jcfpm``, and :yref:`on joint<JCFpmPhys.isOnJoint>` or :yref:`cohesive<JCFpmPhys.isCohesive>` interactions can be vizualized.\n\t``cracks``\n\t\tSaves other data pertaining to the :yref:`rock model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``cracks`` shows locations where cohesive bonds failed during the simulation, with their types (0/1 for tensile/shear breakages), their sizes (0.5*(R1+R2)), and their normal directions. The :yref:`corresponding attribute<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>` has to be activated, and Key attributes have to be consistent.\n\n"))
+ ((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``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n\t``stress``\n\t\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n\t``pericell``\n\t\tSaves the shape of the cell (simulation has to be periodic).\n\t``bstresses``\n\t\tSaves per-particle principal stresses (sigI >= sigII >= sigIII) and associated principal directions (dirI/II/III). Per-particle stress tensor is computed from :yref:`bodyStressTensors<yade.utils.bodyStressTensors>` considering positiv values for compressive states.\n\n.. admonition:: Specific recorders\n\n\tThe following should only be activated in appropriate cases, otherwise crashes can occur due to violation of type presuppositions.\n\n\t``cpm``\n\t\tSaves data pertaining to the :yref:`concrete model<Law2_ScGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmStress`` (stress on particle); ``intr`` is activated automatically by ``cpm``\n\t``wpm``\n\t\tSaves data pertaining to the :yref:`wire particle model<Law2_ScGeom_WirePhys_WirePM>`: ``wpmForceNFactor`` shows the loading factor for the wire, e.g. normal force divided by threshold normal force.\n\t``jcfpm``\n\t\tSaves data pertaining to the :yref:`rock (smooth)-jointed model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``damage`` is defined by :yref:`JCFpmState.tensBreak` + :yref:`JCFpmState.shearBreak`; ``intr`` is activated automatically by ``jcfpm``, and :yref:`on joint<JCFpmPhys.isOnJoint>` or :yref:`cohesive<JCFpmPhys.isCohesive>` interactions can be vizualized.\n\t``cracks``\n\t\tSaves other data pertaining to the :yref:`rock model<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM>`: ``cracks`` shows locations where cohesive bonds failed during the simulation, with their types (0/1 for tensile/shear breakages), their sizes (0.5*(R1+R2)), and their normal directions. The :yref:`corresponding attribute<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>` has to be activated, and Key attributes have to be consistent.\n\n"))
((string,Key,"",,"Necessary if :yref:`recorders<VTKRecorder.recorders>` contains 'cracks'. A string specifying the name of file 'cracks___.txt' that is considered in this case (see :yref:`corresponding attribute<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.Key>`)."))
((int,mask,0,,"If mask defined, only bodies with corresponding groupMask will be exported. If 0, all bodies will be exported.")),
/*ctor*/
=== modified file 'pkg/lbm/HydrodynamicsLawLBM.cpp'
--- pkg/lbm/HydrodynamicsLawLBM.cpp 2014-07-03 07:16:58 +0000
+++ pkg/lbm/HydrodynamicsLawLBM.cpp 2014-07-14 20:08:34 +0000
@@ -26,10 +26,9 @@
#include"HydrodynamicsLawLBM.hpp"
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-#include<boost/filesystem/operations.hpp>
#include<yade/pkg/common/Box.hpp>
#include<yade/pkg/common/Sphere.hpp>
-#include<yade/pkg/dem/CohFrictPhys.hpp>
+#include<yade/pkg/dem/CohesiveFrictionalContactLaw.hpp>
namespace bfs=boost::filesystem;
=== modified file 'py/mathWrap/miniEigen.cpp'
--- py/mathWrap/miniEigen.cpp 2014-07-03 17:24:54 +0000
+++ py/mathWrap/miniEigen.cpp 2014-07-11 17:39:21 +0000
@@ -146,8 +146,8 @@
EIG_WRAP_METH0(Vector2i,Zero); EIG_WRAP_METH0(Vector2i,UnitX); EIG_WRAP_METH0(Vector2i,UnitY); EIG_WRAP_METH0(Vector2i,Ones);
EIG_WRAP_METH0(Quaternionr,Identity);
-#define EIG_OP1(klass,op,sym) TYPEOF((sym klass()).eval()) klass##op(const klass& self){ return (sym self).eval();}
-#define EIG_OP2(klass,op,sym,klass2) TYPEOF((klass() sym klass2()).eval()) klass##op##klass2(const klass& self, const klass2& other){ return (self sym other).eval(); }
+#define EIG_OP1(klass,op,sym) decltype((sym klass()).eval()) klass##op(const klass& self){ return (sym self).eval();}
+#define EIG_OP2(klass,op,sym,klass2) decltype((klass() sym klass2()).eval()) klass##op##klass2(const klass& self, const klass2& other){ return (self sym other).eval(); }
#define EIG_OP2_INPLACE(klass,op,sym,klass2) klass klass##op##klass2(klass& self, const klass2& other){ self sym other; return self; }
=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp 2014-07-03 17:24:54 +0000
+++ py/wrapper/yadeWrapper.cpp 2014-07-14 20:08:34 +0000
@@ -9,7 +9,6 @@
#include<boost/bind.hpp>
#include<boost/lambda/bind.hpp>
#include<boost/thread/thread.hpp>
-#include<boost/filesystem/operations.hpp>
#include<boost/date_time/posix_time/posix_time.hpp>
#include<boost/algorithm/string.hpp>
@@ -683,7 +682,7 @@
void switchToScene(int i){OMEGA.switchToScene(i);}
string sceneToString(){
ostringstream oss;
- yade::ObjectIO::save<TYPEOF(OMEGA.getScene()),boost::archive::binary_oarchive>(oss,"scene",OMEGA.getScene());
+ yade::ObjectIO::save<decltype(OMEGA.getScene()),boost::archive::binary_oarchive>(oss,"scene",OMEGA.getScene());
oss.flush();
return oss.str();
}
