yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11430
[Branch ~yade-pkg/yade/git-trunk] Rev 3426: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Bruno Chareyre (bruno-chareyre)
Jérôme Duriez (jduriez)
Kneib François (francois-kneib)
------------------------------------------------------------
revno: 3426 [merge]
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Wed 2014-05-28 17:12:53 +0200
message:
Merge github.com:yade/trunk into chaoUnsat
added:
extra/
extra/floating_point_utilities_v3/
extra/floating_point_utilities_v3/README
extra/floating_point_utilities_v3/boost/
extra/floating_point_utilities_v3/boost/math/
extra/floating_point_utilities_v3/boost/math/detail/
extra/floating_point_utilities_v3/boost/math/detail/fp_traits.hpp
extra/floating_point_utilities_v3/boost/math/fpclassify.hpp
extra/floating_point_utilities_v3/boost/math/nonfinite_num_facets.hpp
extra/floating_point_utilities_v3/boost/math/signbit.hpp
scripts/checks-and-tests/checks/checkLiquidMigration.py
modified:
CMakeLists.txt
core/main/main.py.in
core/main/yade-batch.in
doc/sphinx/prog.rst
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
pkg/dem/CohFrictMat.hpp
pkg/dem/CohFrictPhys.hpp
pkg/dem/CohesiveFrictionalContactLaw.cpp
pkg/dem/CohesiveFrictionalContactLaw.hpp
pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp
pkg/dem/JointedCohesiveFrictionalPM.cpp
pkg/dem/JointedCohesiveFrictionalPM.hpp
pkg/dem/KinemSimpleShearBox.cpp
pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp
pkg/dem/NormalInelasticityLaw.hpp
pkg/dem/ViscoelasticCapillarPM.cpp
pkg/dem/ViscoelasticCapillarPM.hpp
pkg/pfv/DFNFlow.cpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.ipp
py/config.py.in
py/ymport.py
scripts/checks-and-tests/checks/DEM-PFV-check.py
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'CMakeLists.txt'
--- CMakeLists.txt 2014-05-23 13:37:50 +0000
+++ CMakeLists.txt 2014-05-27 11:24:16 +0000
@@ -88,7 +88,7 @@
# Add possibility to use local boost installation (e.g. -DLocalBoost=1.46.1)
IF ( NOT LocalBoost )
- SET(LocalBoost "1.47.0") # Minimal required Boost version
+ SET(LocalBoost "1.35.0") # Minimal required Boost version
ENDIF ( NOT LocalBoost )
FIND_PACKAGE(Boost ${LocalBoost} COMPONENTS python thread filesystem iostreams regex serialization system date_time REQUIRED)
INCLUDE_DIRECTORIES (${Boost_INCLUDE_DIRS})
@@ -157,7 +157,6 @@
INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS})
LINK_DIRECTORIES( ${VTK_LIBRARY_DIRS} )
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DYADE_VTK")
- SET(features "${features} vtk")
MESSAGE(STATUS "Found VTK")
SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} VTK")
ELSE(VTK_FOUND)
@@ -175,7 +174,6 @@
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DYADE_OPENMP ${OpenMP_CXX_FLAGS}")
MESSAGE(STATUS "Found OpenMP")
SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} OpenMP")
- SET(features "${features} openmp")
ELSE(OPENMP_FOUND)
MESSAGE(STATUS "OpenMP NOT found")
SET(ENABLE_OPENMP OFF)
@@ -195,7 +193,6 @@
INCLUDE_DIRECTORIES(${GLIB2_INCLUDE_DIRS})
MESSAGE(STATUS "Found GTS")
SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} GTS")
- SET(features "${features} gts")
ELSE(GTS_FOUND AND GLIB2_FOUND)
MESSAGE(STATUS "GTS NOT found")
SET(DISABLED_FEATS "${DISABLED_FEATS} GTS")
@@ -214,7 +211,6 @@
IF(QT4_FOUND AND OPENGL_FOUND AND GLUT_FOUND AND GLIB2_FOUND AND QGLVIEWER_FOUND)
SET(GUI_LIBS ${GLUT_LIBRARY} ${OPENGL_LIBRARY} ${QGLVIEWER_LIBRARIES})
SET(GUI_SRC_LIB "lib/opengl/GLUtils.cpp")
- SET(features "${features} qt4 opengl")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DYADE_OPENGL")
INCLUDE_DIRECTORIES(${GLIB2_INCLUDE_DIRS})
INCLUDE_DIRECTORIES(${QT_INCLUDES})
@@ -255,7 +251,7 @@
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")
ENDIF(ENABLE_PFVFLOW)
@@ -330,6 +326,16 @@
INCLUDE_DIRECTORIES(${CMAKE_BINARY_DIR})
#===========================================================
+# floating_point_utilities_v3 are already in Boost included
+# Use embedded copy only if Boost older than 1.47.0
+
+IF((Boost_MINOR_VERSION LESS 47) AND (Boost_MAJOR_VERSION EQUAL 1))
+ MESSAGE(STATUS "Boost version is less than 1.47, using embedded version of floating_point_utilities_v3")
+ MESSAGE(STATUS "Consider updating boost or system, as this embedded library will be removed soon")
+ INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/extra/floating_point_utilities_v3)
+ENDIF((Boost_MINOR_VERSION LESS 47) AND (Boost_MAJOR_VERSION EQUAL 1))
+
+#===========================================================
IF(ENABLE_LBMFLOW)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DLBM_ENGINE")
SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} LBMFLOW")
=== modified file 'core/main/main.py.in'
--- core/main/main.py.in 2014-02-03 16:10:20 +0000
+++ core/main/main.py.in 2014-05-27 11:24:16 +0000
@@ -18,7 +18,9 @@
# get yade path (allow YADE_PREFIX to override)
prefix,suffix='${runtimePREFIX}' if not os.environ.has_key('YADE_PREFIX') else os.environ['YADE_PREFIX'],'${SUFFIX}'
# duplicate some items from yade.config here, so that we can increase verbosity when the c++ part is booting
-features,version,debugbuild='${features}'.split(' '),'${realVersion}',' ${debugbuild}'
+features,version,debugbuild='${CONFIGURED_FEATS}'.split(' '),'${realVersion}',' ${debugbuild}'
+
+if (features[0]==''): features=features[1:]
libPATH='${LIBRARY_OUTPUT_PATH}'
if (libPATH[1:] == '{LIBRARY_OUTPUT_PATH}'): libPATH='lib'
@@ -101,7 +103,7 @@
print 20*'*'+' SOME TESTS FAILED '+20*'*'
sys.exit(1)
-if not 'openmp' in features and (opts.cores or (opts.threads and opts.threads>1)):
+if not 'OpenMP' in features and (opts.cores or (opts.threads and opts.threads>1)):
print 'WARNING: compiled without OpenMP, -j/--threads/--cores have no effect.'
# OpenMP env variables must be set before loading yade libs ("import yade" below)
@@ -240,7 +242,7 @@
## run userSession in a way corresponding to the features we use:
gui=None
yade.runtime.hasDisplay=False # this is the default initialized in the module, anyway
-if 'qt4' in features: gui='qt4'
+if 'GUI' in features: gui='qt4'
if opts.nogui: gui=None
if gui:
import Xlib.display
=== modified file 'core/main/yade-batch.in'
--- core/main/yade-batch.in 2013-12-20 15:10:54 +0000
+++ core/main/yade-batch.in 2014-05-27 11:24:16 +0000
@@ -343,8 +343,9 @@
if os.environ.has_key("OMP_NUM_THREADS"): return min(int(os.environ['OMP_NUM_THREADS']),nCpu)
return nCpu
numCores=getNumCores()
-maxOmpThreads=numCores if 'openmp' in yade.config.features else 1
-features,version='${features}'.split(','),'${realVersion}'
+maxOmpThreads=numCores if 'OpenMP' in yade.config.features else 1
+features,version='${CONFIGURED_FEATS}'.split(','),'${realVersion}'
+if (features[0]==''): features=features[1:]
prog = os.path.basename(sys.argv[0])
parser=argparse.ArgumentParser(usage='%s [options] [ TABLE [SIMULATION.py] | SIMULATION.py[/nCores] [...] ]'%prog,description='%s runs yade simulation multiple times with different parameters.\n\nSee https://yade-dem.org/sphinx/user.html#batch-queuing-and-execution-yade-batch for details.\n\nBatch can be specified either with parameter table TABLE (must not end in .py), which is either followed by exactly one SIMULATION.py (must end in .py), or contains !SCRIPT column specifying the simulation to be run. The second option is to specify multiple scripts, which can optionally have /nCores suffix to specify number of cores for that particular simulation (corresponds to !THREADS column in the parameter table), e.g. sim.py/3.'%prog)
@@ -473,7 +474,7 @@
nCores=maxJobs
else:
logging.warning('WARNING: job #%d will use %d cores but only %d are available'%(i,nCores,maxJobs))
- if 'openmp' not in yade.config.features and nCores>1:
+ if 'OpenMP' not in yade.config.features and nCores>1:
logging.warning('Job #%d will be uselessly run with %d threads (compiled without OpenMP support).'%(i,nCores))
executables.add(jobExecutable)
# compose command-line: build the hyper-linked variant, then strip HTML tags (ugly, but ensures consistency)
=== modified file 'doc/sphinx/prog.rst'
--- doc/sphinx/prog.rst 2013-09-12 06:56:25 +0000
+++ doc/sphinx/prog.rst 2014-05-27 11:24:16 +0000
@@ -62,7 +62,7 @@
YADE_REQUIRE_FEATURE(vtk);
YADE_REQUIRE_FEATURE(gts);
-This file will be compiled only if *both* ``vtk`` and ``gts`` features are enabled.
+This file will be compiled only if *both* ``VTK`` and ``GTS`` features are enabled.
Depending on current feature set, only selection of plugins will be compiled.
It is possible to disable compilation of a file by requiring any non-existent feature, such as::
=== added directory 'extra'
=== added directory 'extra/floating_point_utilities_v3'
=== added file 'extra/floating_point_utilities_v3/README'
--- extra/floating_point_utilities_v3/README 1970-01-01 00:00:00 +0000
+++ extra/floating_point_utilities_v3/README 2014-05-26 08:10:59 +0000
@@ -0,0 +1,7 @@
+Graceful handling of nan and inf in iostream
+
+Proposed for boost::math
+
+Downloaded from: http://www.boostpro.com/vault/index.php?action=downloadfile&filename=floating_point_utilities_v3.zip&directory=Math%20-%20Numerics&;
+
+stripped of documentation
=== added directory 'extra/floating_point_utilities_v3/boost'
=== added directory 'extra/floating_point_utilities_v3/boost/math'
=== added directory 'extra/floating_point_utilities_v3/boost/math/detail'
=== added file 'extra/floating_point_utilities_v3/boost/math/detail/fp_traits.hpp'
--- extra/floating_point_utilities_v3/boost/math/detail/fp_traits.hpp 1970-01-01 00:00:00 +0000
+++ extra/floating_point_utilities_v3/boost/math/detail/fp_traits.hpp 2014-05-26 08:10:59 +0000
@@ -0,0 +1,576 @@
+// fp_traits.hpp
+
+#ifndef BOOST_MATH_FP_TRAITS_HPP
+#define BOOST_MATH_FP_TRAITS_HPP
+
+// Copyright (c) 2006 Johan Rade
+
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#if defined(__vms) && defined(__DECCXX) && !__IEEE_FLOAT
+# error The VAX floating point mode on VMS is not supported.
+#endif
+
+#include <cstring>
+
+#include <boost/assert.hpp>
+#include <boost/cstdint.hpp>
+#include <boost/detail/endian.hpp>
+#include <boost/static_assert.hpp>
+#include <boost/type_traits/is_floating_point.hpp>
+
+//------------------------------------------------------------------------------
+
+namespace boost {
+namespace math {
+namespace detail {
+
+//------------------------------------------------------------------------------
+
+/*
+Most processors support three different floating point precisions:
+single precision (32 bits), double precision (64 bits)
+and extended double precision (>64 bits)
+
+Note that the C++ type long double can be implemented
+both as double precision and extended double precision.
+*/
+
+struct single_precision_tag {};
+struct double_precision_tag {};
+struct extended_double_precision_tag {};
+
+//------------------------------------------------------------------------------
+
+/*
+template<class T, class U> struct fp_traits_impl;
+
+ This is traits class that describes the binary structure of floating
+ point numbers of C++ type T and precision U
+
+Requirements:
+
+ T = float, double or long double
+ U = single_precision_tag, double_precision_tag
+ or extended_double_precision_tag
+
+Typedef members:
+
+ bits -- the target type when copying the leading bytes of a floating
+ point number. It is a typedef for uint32_t or uint64_t.
+
+ coverage -- tells us whether all bytes are copied or not.
+ It is a typedef for all_bits or not_all_bits.
+
+Static data members:
+
+ sign, exponent, flag, mantissa -- bit masks that give the meaning of the bits
+ in the leading bytes.
+
+Static function members:
+
+ init() -- initializes the static data members, if needed.
+ (Is a no-op in the specialized versions of the template.)
+
+ get_bits(), set_bits() -- provide access to the leading bytes.
+*/
+
+struct all_bits {};
+struct not_all_bits {};
+
+// Generic version -------------------------------------------------------------
+
+// The generic version uses run time initialization to determine the floating
+// point format. It is capable of handling most formats,
+// but not the Motorola 68K extended double precision format.
+
+// Currently the generic version is used only for extended double precision
+// on Itanium. In all other cases there are specializations of the template
+// that use compile time initialization.
+
+template<class T> struct uint32_t_coverage
+{
+ typedef not_all_bits type;
+};
+
+template<> struct uint32_t_coverage<single_precision_tag>
+{
+ typedef all_bits type;
+};
+
+template<class T, class U> struct fp_traits_impl
+{
+ typedef uint32_t bits;
+ typedef BOOST_DEDUCED_TYPENAME uint32_t_coverage<U>::type coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ static uint32_t exponent;
+ static uint32_t flag;
+ static uint32_t mantissa;
+
+ static void init()
+ {
+ if(is_init_) return;
+ do_init_();
+ is_init_ = true;
+ }
+
+ static void get_bits(T x, uint32_t& a)
+ {
+ memcpy(&a, reinterpret_cast<const unsigned char*>(&x) + offset_, 4);
+ }
+
+ static void set_bits(T& x, uint32_t a)
+ {
+ memcpy(reinterpret_cast<unsigned char*>(&x) + offset_, &a, 4);
+ }
+
+private:
+ static size_t offset_;
+ static bool is_init_;
+ static void do_init_();
+};
+
+//..............................................................................
+
+template<class T, class U> uint32_t fp_traits_impl<T,U>::exponent;
+template<class T, class U> uint32_t fp_traits_impl<T,U>::flag;
+template<class T, class U> uint32_t fp_traits_impl<T,U>::mantissa;
+template<class T, class U> size_t fp_traits_impl<T,U>::offset_;
+template<class T, class U> bool fp_traits_impl<T,U>::is_init_;
+
+// In a single-threaded program, do_init will be called exactly once.
+// In a multi-threaded program, do_init may be called simultaneously
+// by more then one thread. That should not be a problem.
+
+//..............................................................................
+
+template<class T, class U> void fp_traits_impl<T,U>::do_init_()
+{
+ T x = static_cast<T>(3) / static_cast<T>(4);
+ // sign bit = 0
+ // exponent: first and last bit = 0, all other bits = 1
+ // flag bit (if present) = 1
+ // mantissa: first bit = 1, all other bits = 0
+
+ uint32_t a;
+
+ for(size_t k = 0; k <= sizeof(T) - 4; ++k) {
+
+ memcpy(&a, reinterpret_cast<unsigned char*>(&x) + k, 4);
+
+ switch(a) {
+
+ case 0x3f400000: // IEEE single precision format
+
+ offset_ = k;
+ exponent = 0x7f800000;
+ flag = 0x00000000;
+ mantissa = 0x007fffff;
+ return;
+
+ case 0x3fe80000: // IEEE double precision format
+ // and PowerPC extended double precision format
+ offset_ = k;
+ exponent = 0x7ff00000;
+ flag = 0x00000000;
+ mantissa = 0x000fffff;
+ return;
+
+ case 0x3ffe0000: // Motorola extended double precision format
+
+ // Must not get here. Must be handled by specialization.
+ // To get accurate cutoff between normals and subnormals
+ // we must use the flag bit that is in the 5th byte.
+ // Otherwise this cutoff will be off by a factor 2.
+ // If we do get here, then we have failed to detect the Motorola
+ // processor at compile time.
+
+ BOOST_ASSERT(false);
+ return;
+
+ case 0x3ffe8000: // IEEE extended double precision format
+ // with 15 exponent bits
+ offset_ = k;
+ exponent = 0x7fff0000;
+ flag = 0x00000000;
+ mantissa = 0x0000ffff;
+ return;
+
+ case 0x3ffec000: // Intel extended double precision format
+
+ offset_ = k;
+ exponent = 0x7fff0000;
+ flag = 0x00008000;
+ mantissa = 0x00007fff;
+ return;
+
+ default:
+ continue;
+ }
+ }
+
+ BOOST_ASSERT(false);
+
+ // Unknown format.
+}
+
+
+// float (32 bits) -------------------------------------------------------------
+
+template<> struct fp_traits_impl<float, single_precision_tag>
+{
+ typedef uint32_t bits;
+ typedef all_bits coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ BOOST_STATIC_CONSTANT(uint32_t, exponent = 0x7f800000);
+ BOOST_STATIC_CONSTANT(uint32_t, flag = 0x00000000);
+ BOOST_STATIC_CONSTANT(uint32_t, mantissa = 0x007fffff);
+
+ static void init() {}
+ static void get_bits(float x, uint32_t& a) { memcpy(&a, &x, 4); }
+ static void set_bits(float& x, uint32_t a) { memcpy(&x, &a, 4); }
+};
+
+
+// double (64 bits) ------------------------------------------------------------
+
+#if defined(BOOST_NO_INT64_T) || defined(BOOST_NO_INCLASS_MEMBER_INITIALIZATION)
+
+template<> struct fp_traits_impl<double, double_precision_tag>
+{
+ typedef uint32_t bits;
+ typedef not_all_bits coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ BOOST_STATIC_CONSTANT(uint32_t, exponent = 0x7ff00000);
+ BOOST_STATIC_CONSTANT(uint32_t, flag = 0);
+ BOOST_STATIC_CONSTANT(uint32_t, mantissa = 0x000fffff);
+
+ static void init() {}
+
+ static void get_bits(double x, uint32_t& a)
+ {
+ memcpy(&a, reinterpret_cast<const unsigned char*>(&x) + offset_, 4);
+ }
+
+ static void set_bits(double& x, uint32_t a)
+ {
+ memcpy(reinterpret_cast<unsigned char*>(&x) + offset_, &a, 4);
+ }
+
+private:
+
+#if defined(BOOST_BIG_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 0);
+#elif defined(BOOST_LITTLE_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 4);
+#else
+ BOOST_STATIC_ASSERT(false);
+#endif
+};
+
+//..............................................................................
+
+#else
+
+template<> struct fp_traits_impl<double, double_precision_tag>
+{
+ typedef uint64_t bits;
+ typedef all_bits coverage;
+
+ static const uint64_t sign = (uint64_t)0x80000000 << 32;
+ static const uint64_t exponent = (uint64_t)0x7ff00000 << 32;
+ static const uint64_t flag = 0;
+ static const uint64_t mantissa
+ = ((uint64_t)0x000fffff << 32) + (uint64_t)0xffffffff;
+
+ static void init() {}
+ static void get_bits(double x, uint64_t& a) { memcpy(&a, &x, 8); }
+ static void set_bits(double& x, uint64_t a) { memcpy(&x, &a, 8); }
+};
+
+#endif
+
+
+// long double (64 bits) -------------------------------------------------------
+
+#if defined(BOOST_NO_INT64_T) || defined(BOOST_NO_INCLASS_MEMBER_INITIALIZATION)
+
+template<> struct fp_traits_impl<long double, double_precision_tag>
+{
+ typedef uint32_t bits;
+ typedef not_all_bits coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ BOOST_STATIC_CONSTANT(uint32_t, exponent = 0x7ff00000);
+ BOOST_STATIC_CONSTANT(uint32_t, flag = 0);
+ BOOST_STATIC_CONSTANT(uint32_t, mantissa = 0x000fffff);
+
+ static void init() {}
+
+ static void get_bits(long double x, uint32_t& a)
+ {
+ memcpy(&a, reinterpret_cast<const unsigned char*>(&x) + offset_, 4);
+ }
+
+ static void set_bits(long double& x, uint32_t a)
+ {
+ memcpy(reinterpret_cast<unsigned char*>(&x) + offset_, &a, 4);
+ }
+
+private:
+
+#if defined(BOOST_BIG_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 0);
+#elif defined(BOOST_LITTLE_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 4);
+#else
+ BOOST_STATIC_ASSERT(false);
+#endif
+};
+
+//..............................................................................
+
+#else
+
+template<> struct fp_traits_impl<long double, double_precision_tag>
+{
+ typedef uint64_t bits;
+ typedef all_bits coverage;
+
+ static const uint64_t sign = (uint64_t)0x80000000 << 32;
+ static const uint64_t exponent = (uint64_t)0x7ff00000 << 32;
+ static const uint64_t flag = 0;
+ static const uint64_t mantissa
+ = ((uint64_t)0x000fffff << 32) + (uint64_t)0xffffffff;
+
+ static void init() {}
+ static void get_bits(long double x, uint64_t& a) { memcpy(&a, &x, 8); }
+ static void set_bits(long double& x, uint64_t a) { memcpy(&x, &a, 8); }
+};
+
+#endif
+
+
+// long double (>64 bits), x86 and x64 -----------------------------------------
+
+#if defined(__i386) || defined(__i386__) || defined(_M_IX86) \
+ || defined(__amd64) || defined(__amd64__) || defined(_M_AMD64) \
+ || defined(__x86_64) || defined(__x86_64__) || defined(_M_X64)
+
+// Intel extended double precision format (80 bits)
+
+template<> struct fp_traits_impl<long double, extended_double_precision_tag>
+{
+ typedef uint32_t bits;
+ typedef not_all_bits coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ BOOST_STATIC_CONSTANT(uint32_t, exponent = 0x7fff0000);
+ BOOST_STATIC_CONSTANT(uint32_t, flag = 0x00008000);
+ BOOST_STATIC_CONSTANT(uint32_t, mantissa = 0x00007fff);
+
+ static void init() {}
+
+ static void get_bits(long double x, uint32_t& a)
+ {
+ memcpy(&a, reinterpret_cast<const unsigned char*>(&x) + 6, 4);
+ }
+
+ static void set_bits(long double& x, uint32_t a)
+ {
+ memcpy(reinterpret_cast<unsigned char*>(&x) + 6, &a, 4);
+ }
+};
+
+
+// long double (>64 bits), Itanium ---------------------------------------------
+
+#elif defined(__ia64) || defined(__ia64__) || defined(_M_IA64)
+
+// The floating point format is unknown at compile time
+// No template specialization is provided.
+// The generic definition is used.
+
+// The Itanium supports both
+// the Intel extended double precision format (80 bits) and
+// the IEEE extended double precision format with 15 exponent bits (128 bits).
+
+
+// long double (>64 bits), PowerPC ---------------------------------------------
+
+#elif defined(__powerpc) || defined(__powerpc__) || defined(__POWERPC__) \
+ || defined(__ppc) || defined(__ppc__) || defined(__PPC__)
+
+// PowerPC extended double precision format (128 bits)
+
+template<> struct fp_traits_impl<long double, extended_double_precision_tag>
+{
+ typedef uint32_t bits;
+ typedef not_all_bits coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ BOOST_STATIC_CONSTANT(uint32_t, exponent = 0x7ff00000);
+ BOOST_STATIC_CONSTANT(uint32_t, flag = 0x00000000);
+ BOOST_STATIC_CONSTANT(uint32_t, mantissa = 0x000fffff);
+
+ static void init() {}
+
+ static void get_bits(long double x, uint32_t& a)
+ {
+ memcpy(&a, reinterpret_cast<const unsigned char*>(&x) + offset_, 4);
+ }
+
+ static void set_bits(long double& x, uint32_t a)
+ {
+ memcpy(reinterpret_cast<unsigned char*>(&x) + offset_, &a, 4);
+ }
+
+private:
+
+#if defined(BOOST_BIG_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 0);
+#elif defined(BOOST_LITTLE_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 12);
+#else
+ BOOST_STATIC_ASSERT(false);
+#endif
+};
+
+
+// long double (>64 bits), Motorola 68K ----------------------------------------
+
+#elif defined(__m68k) || defined(__m68k__) \
+ || defined(__mc68000) || defined(__mc68000__) \
+
+// Motorola extended double precision format (96 bits)
+
+// It is the same format as the Intel extended double precision format,
+// except that 1) it is big-endian, 2) the 3rd and 4th byte are padding, and
+// 3) the flag bit is not set for infinity
+
+template<> struct fp_traits_impl<long double, extended_double_precision_tag>
+{
+ typedef uint32_t bits;
+ typedef not_all_bits coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ BOOST_STATIC_CONSTANT(uint32_t, exponent = 0x7fff0000);
+ BOOST_STATIC_CONSTANT(uint32_t, flag = 0x00008000);
+ BOOST_STATIC_CONSTANT(uint32_t, mantissa = 0x00007fff);
+
+ static void init() {}
+
+ // copy 1st, 2nd, 5th and 6th byte. 3rd and 4th byte are padding.
+
+ static void get_bits(long double x, uint32_t& a)
+ {
+ memcpy(&a, &x, 2);
+ memcpy(reinterpret_cast<unsigned char*>(&a) + 2,
+ reinterpret_cast<const unsigned char*>(&x) + 4, 2);
+ }
+
+ static void set_bits(long double& x, uint32_t a)
+ {
+ memcpy(&x, &a, 2);
+ memcpy(reinterpret_cast<unsigned char*>(&x) + 4,
+ reinterpret_cast<const unsigned char*>(&a) + 2, 2);
+ }
+};
+
+
+// long double (>64 bits), All other processors --------------------------------
+
+#else
+
+// IEEE extended double precision format with 15 exponent bits (128 bits)
+
+template<> struct fp_traits_impl<long double, extended_double_precision_tag>
+{
+ typedef uint32_t bits;
+ typedef not_all_bits coverage;
+
+ BOOST_STATIC_CONSTANT(uint32_t, sign = 0x80000000);
+ BOOST_STATIC_CONSTANT(uint32_t, exponent = 0x7fff0000);
+ BOOST_STATIC_CONSTANT(uint32_t, flag = 0x00000000);
+ BOOST_STATIC_CONSTANT(uint32_t, mantissa = 0x0000ffff);
+
+ static void init() {}
+
+ static void get_bits(long double x, uint32_t& a)
+ {
+ memcpy(&a, reinterpret_cast<const unsigned char*>(&x) + offset_, 4);
+ }
+
+ static void set_bits(long double& x, uint32_t a)
+ {
+ memcpy(reinterpret_cast<unsigned char*>(&x) + offset_, &a, 4);
+ }
+
+private:
+
+#if defined(BOOST_BIG_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 0);
+#elif defined(BOOST_LITTLE_ENDIAN)
+ BOOST_STATIC_CONSTANT(int, offset_ = 12);
+#else
+ BOOST_STATIC_ASSERT(false);
+#endif
+};
+
+#endif
+
+
+//------------------------------------------------------------------------------
+
+// size_to_precision is a type switch for converting a C++ floating point type
+// to the corresponding precision type.
+
+template<int n> struct size_to_precision;
+
+template<> struct size_to_precision<4>
+{
+ typedef single_precision_tag type;
+};
+
+template<> struct size_to_precision<8>
+{
+ typedef double_precision_tag type;
+};
+
+template<> struct size_to_precision<10>
+{
+ typedef extended_double_precision_tag type;
+};
+
+template<> struct size_to_precision<12>
+{
+ typedef extended_double_precision_tag type;
+};
+
+template<> struct size_to_precision<16>
+{
+ typedef extended_double_precision_tag type;
+};
+
+// fp_traits is a type switch that selects the right fp_traits_impl
+
+template<class T> struct fp_traits
+{
+ BOOST_STATIC_ASSERT(boost::is_floating_point<T>::value);
+ typedef BOOST_DEDUCED_TYPENAME size_to_precision<sizeof(T)>::type precision;
+ typedef fp_traits_impl<T, precision> type;
+};
+
+
+//------------------------------------------------------------------------------
+
+} // namespace detail
+} // namespace math
+} // namespace boost
+
+#endif
=== added file 'extra/floating_point_utilities_v3/boost/math/fpclassify.hpp'
--- extra/floating_point_utilities_v3/boost/math/fpclassify.hpp 1970-01-01 00:00:00 +0000
+++ extra/floating_point_utilities_v3/boost/math/fpclassify.hpp 2014-05-26 08:10:59 +0000
@@ -0,0 +1,229 @@
+// fpclassify.hpp
+
+#ifndef BOOST_MATH_FPCLASSIFY_HPP
+#define BOOST_MATH_FPCLASSIFY_HPP
+
+// Copyright (c) 2006 Johan Rade
+
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+/*
+The following algorithm is used:
+
+ If all exponent bits, the flag bit (if there is one),
+ and all mantissa bits are 0, then the number is zero.
+
+ If all exponent bits and the flag bit (if there is one) are 0,
+ and at least one mantissa bit is 1, then the number is subnormal.
+
+ If all exponent bits are 1 and all mantissa bits are 0,
+ then the number is infinity.
+
+ If all exponent bits are 1 and at least one mantissa bit is 1,
+ then the number is a not-a-number.
+
+ Otherwise the number is normal.
+
+(Note that the binary representation of infinity
+has flag bit 0 for Motorola 68K extended double precision,
+and flag bit 1 for Intel extended double precision.)
+
+To get the bits, the four or eight most significant bytes are copied
+into an uint32_t or uint64_t and bit masks are applied.
+This covers all the exponent bits and the flag bit (if there is one),
+but not always all the mantissa bits.
+Some of the functions below have two implementations,
+depending on whether all the mantissa bits are copied or not.
+*/
+
+#include <cmath>
+
+#ifndef FP_INFINITE
+# define FP_INFINITE 0
+# define FP_NAN 1
+# define FP_NORMAL 2
+# define FP_SUBNORMAL 3
+# define FP_ZERO 4
+#endif
+
+#include "detail/fp_traits.hpp"
+
+namespace boost {
+namespace math {
+
+//------------------------------------------------------------------------------
+
+template<class T> bool (isfinite)(T x)
+{
+ typedef BOOST_DEDUCED_TYPENAME detail::fp_traits<T>::type traits;
+ traits::init();
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::exponent;
+ return a != traits::exponent;
+}
+
+//------------------------------------------------------------------------------
+
+template<class T> bool (isnormal)(T x)
+{
+ typedef BOOST_DEDUCED_TYPENAME detail::fp_traits<T>::type traits;
+ traits::init();
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::exponent | traits::flag;
+ return (a != 0) && (a < traits::exponent);
+}
+
+//------------------------------------------------------------------------------
+
+namespace detail {
+
+ template<class T> bool isinf_impl(T x, all_bits)
+ {
+ typedef BOOST_DEDUCED_TYPENAME fp_traits<T>::type traits;
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::exponent | traits::mantissa;
+ return a == traits::exponent;
+ }
+
+ template<class T> bool isinf_impl(T x, not_all_bits)
+ {
+ typedef BOOST_DEDUCED_TYPENAME fp_traits<T>::type traits;
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::exponent | traits::mantissa;
+ if(a != traits::exponent)
+ return false;
+
+ traits::set_bits(x,0);
+ return x == 0;
+ }
+
+} // namespace detail
+
+template<class T> bool (isinf)(T x)
+{
+ typedef BOOST_DEDUCED_TYPENAME detail::fp_traits<T>::type traits;
+ traits::init();
+ return detail::isinf_impl(x, BOOST_DEDUCED_TYPENAME traits::coverage());
+}
+
+//------------------------------------------------------------------------------
+
+namespace detail {
+
+ template<class T> bool isnan_impl(T x, all_bits)
+ {
+ typedef BOOST_DEDUCED_TYPENAME fp_traits<T>::type traits;
+ traits::init();
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::exponent | traits::mantissa;
+ return a > traits::exponent;
+ }
+
+ template<class T> bool isnan_impl(T x, not_all_bits)
+ {
+ typedef BOOST_DEDUCED_TYPENAME fp_traits<T>::type traits;
+ traits::init();
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+
+ a &= traits::exponent | traits::mantissa;
+ if(a < traits::exponent)
+ return false;
+
+ a &= traits::mantissa;
+ traits::set_bits(x,a);
+ return x != 0;
+ }
+
+} // namespace detail
+
+template<class T> bool (isnan)(T x)
+{
+ typedef BOOST_DEDUCED_TYPENAME detail::fp_traits<T>::type traits;
+ traits::init();
+ return detail::isnan_impl(x, BOOST_DEDUCED_TYPENAME traits::coverage());
+}
+
+//------------------------------------------------------------------------------
+
+namespace detail {
+
+ template<class T> int fpclassify_impl(T x, all_bits)
+ {
+ typedef BOOST_DEDUCED_TYPENAME fp_traits<T>::type traits;
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::exponent | traits::flag | traits::mantissa;
+
+ if(a <= traits::mantissa) {
+ if(a == 0)
+ return FP_ZERO;
+ else
+ return FP_SUBNORMAL;
+ }
+
+ if(a < traits::exponent)
+ return FP_NORMAL;
+
+ a &= traits::mantissa;
+ if(a == 0)
+ return FP_INFINITE;
+
+ return FP_NAN;
+ }
+
+ template<class T> int fpclassify_impl(T x, not_all_bits)
+ {
+ typedef BOOST_DEDUCED_TYPENAME fp_traits<T>::type traits;
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::exponent | traits::flag | traits::mantissa;
+
+ if(a <= traits::mantissa) {
+ if(x == 0)
+ return FP_ZERO;
+ else
+ return FP_SUBNORMAL;
+ }
+
+ if(a < traits::exponent)
+ return FP_NORMAL;
+
+ a &= traits::mantissa;
+ traits::set_bits(x,a);
+ if(x == 0)
+ return FP_INFINITE;
+
+ return FP_NAN;
+ }
+
+} // namespace detail
+
+template<class T> int (fpclassify)(T x)
+{
+ typedef BOOST_DEDUCED_TYPENAME detail::fp_traits<T>::type traits;
+ traits::init();
+ return detail::fpclassify_impl(x, BOOST_DEDUCED_TYPENAME traits::coverage());
+}
+
+//------------------------------------------------------------------------------
+
+} // namespace math
+} // namespace boost
+
+#endif
=== added file 'extra/floating_point_utilities_v3/boost/math/nonfinite_num_facets.hpp'
--- extra/floating_point_utilities_v3/boost/math/nonfinite_num_facets.hpp 1970-01-01 00:00:00 +0000
+++ extra/floating_point_utilities_v3/boost/math/nonfinite_num_facets.hpp 2014-05-26 08:10:59 +0000
@@ -0,0 +1,475 @@
+#ifndef BOOST_MATH_NONFINITE_NUM_FACETS_HPP
+#define BOOST_MATH_NONFINITE_NUM_FACETS_HPP
+
+// Copyright (c) 2006 Johan Rade
+
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <cstring>
+#include <ios>
+#include <limits>
+#include <locale>
+#include "fpclassify.hpp"
+#include "signbit.hpp"
+
+#ifdef _MSC_VER
+# pragma warning(push)
+# pragma warning(disable : 4127 4511 4512 4706)
+#endif
+
+namespace boost {
+namespace math {
+
+
+// flags -----------------------------------------------------------------------
+
+const int legacy = 0x1;
+const int signed_zero = 0x2;
+const int trap_infinity = 0x4;
+const int trap_nan = 0x8;
+
+
+// class nonfinite_num_put -----------------------------------------------------
+
+template<
+ class CharType,
+ class OutputIterator = std::ostreambuf_iterator<CharType>
+>
+class nonfinite_num_put : public std::num_put<CharType, OutputIterator> {
+public:
+ explicit nonfinite_num_put(int flags = 0) : flags_(flags) {}
+
+protected:
+ virtual OutputIterator do_put(
+ OutputIterator it, std::ios_base& iosb,
+ CharType fill, double val) const
+ {
+ put_and_reset_width(it, iosb, fill, val);
+ return it;
+ }
+
+ virtual OutputIterator do_put(
+ OutputIterator it, std::ios_base& iosb,
+ CharType fill, long double val) const
+ {
+ put_and_reset_width(it, iosb, fill, val);
+ return it;
+ }
+
+private:
+ template<class ValType> void put_and_reset_width(
+ OutputIterator& it, std::ios_base& iosb,
+ CharType fill, ValType val) const
+ {
+ put_impl(it, iosb, fill, val);
+ iosb.width(0);
+ }
+
+ template<class ValType> void put_impl(
+ OutputIterator& it, std::ios_base& iosb,
+ CharType fill, ValType val) const
+ {
+ switch((boost::math::fpclassify)(val)) {
+
+ case FP_INFINITE:
+ if(flags_ & trap_infinity)
+ throw std::ios_base::failure("Infinity");
+ else if((boost::math::signbit)(val))
+ put_num_and_fill(it, iosb, "-", "inf", fill);
+ else if(iosb.flags() & std::ios_base::showpos)
+ put_num_and_fill(it, iosb, "+", "inf", fill);
+ else
+ put_num_and_fill(it, iosb, "", "inf", fill);
+ break;
+
+ case FP_NAN:
+ if(flags_ & trap_nan)
+ throw std::ios_base::failure("NaN");
+ else if((boost::math::signbit)(val))
+ put_num_and_fill(it, iosb, "-", "nan", fill);
+ else if(iosb.flags() & std::ios_base::showpos)
+ put_num_and_fill(it, iosb, "+", "nan", fill);
+ else
+ put_num_and_fill(it, iosb, "", "nan", fill);
+ break;
+
+ case FP_ZERO:
+ if(flags_ & signed_zero) {
+ if((boost::math::signbit)(val))
+ put_num_and_fill(it, iosb, "-", "0", fill);
+ else if(iosb.flags() & std::ios_base::showpos)
+ put_num_and_fill(it, iosb, "+", "0", fill);
+ else
+ put_num_and_fill(it, iosb, "", "0", fill);
+ }
+ else
+ put_num_and_fill(it, iosb, "", "0", fill);
+ break;
+
+ default:
+ it = std::num_put<CharType, OutputIterator>::do_put(
+ it, iosb, fill, val);
+ break;
+ }
+ }
+
+ void put_num_and_fill(
+ OutputIterator& it, std::ios_base& iosb, const char* prefix,
+ const char* body, CharType fill) const
+ {
+ int width = (int)strlen(prefix) + (int)strlen(body);
+ std::ios_base::fmtflags adjust
+ = iosb.flags() & std::ios_base::adjustfield;
+ const std::ctype<CharType>& ct
+ = std::use_facet<std::ctype<CharType> >(iosb.getloc());
+
+ if(adjust != std::ios_base::internal && adjust != std::ios_base::left)
+ put_fill(it, iosb, fill, width);
+
+ while(*prefix)
+ *it = ct.widen(*(prefix++));
+
+ if(adjust == std::ios_base::internal)
+ put_fill(it, iosb, fill, width);
+
+ if(iosb.flags() & std::ios_base::uppercase) {
+ while(*body)
+ *it = ct.toupper(ct.widen(*(body++)));
+ }
+ else {
+ while(*body)
+ *it = ct.widen(*(body++));
+ }
+
+ if(adjust == std::ios_base::left)
+ put_fill(it, iosb, fill, width);
+ }
+
+ void put_fill(
+ OutputIterator& it, std::ios_base& iosb,
+ CharType fill, int width) const
+ {
+ for(int i = iosb.width() - width; i > 0; --i)
+ *it = fill;
+ }
+
+private:
+ const int flags_;
+};
+
+
+// class nonfinite_num_get ------------------------------------------------------
+
+template<
+ class CharType,
+ class InputIterator = std::istreambuf_iterator<CharType>
+>
+class nonfinite_num_get : public std::num_get<CharType, InputIterator> {
+public:
+ explicit nonfinite_num_get(int flags = 0) : flags_(flags) {}
+
+protected:
+ virtual InputIterator do_get(
+ InputIterator it, InputIterator end, std::ios_base& iosb,
+ std::ios_base::iostate& state, float& val) const
+ {
+ get_and_check_eof(it, end, iosb, state, val);
+ return it;
+ }
+
+ virtual InputIterator do_get(
+ InputIterator it, InputIterator end, std::ios_base& iosb,
+ std::ios_base::iostate& state, double& val) const
+ {
+ get_and_check_eof(it, end, iosb, state, val);
+ return it;
+ }
+
+ virtual InputIterator do_get(
+ InputIterator it, InputIterator end, std::ios_base& iosb,
+ std::ios_base::iostate& state, long double& val) const
+ {
+ get_and_check_eof(it, end, iosb, state, val);
+ return it;
+ }
+
+//..............................................................................
+
+private:
+ template<class ValType> static ValType positive_nan()
+ {
+ // on some platforms quiet_NaN() is negative
+ return (boost::math::copysign)(
+ std::numeric_limits<ValType>::quiet_NaN(), 1);
+ }
+
+ template<class ValType> void get_and_check_eof(
+ InputIterator& it, InputIterator end, std::ios_base& iosb,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ get_signed(it, end, iosb, state, val);
+ if(it == end)
+ state |= std::ios_base::eofbit;
+ }
+
+ template<class ValType> void get_signed(
+ InputIterator& it, InputIterator end, std::ios_base& iosb,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ const std::ctype<CharType>& ct
+ = std::use_facet<std::ctype<CharType> >(iosb.getloc());
+
+ char c = peek_char(it, end, ct);
+
+ bool negative = (c == '-');
+
+ if(negative || c == '+') {
+ ++it;
+ c = peek_char(it, end, ct);
+ if(c == '-' || c == '+') {
+ // without this check, "++5" etc would be accepted
+ state |= std::ios_base::failbit;
+ return;
+ }
+ }
+
+ get_unsigned(it, end, iosb, ct, state, val);
+
+ if(negative)
+ val = (boost::math::changesign)(val);
+ }
+
+ template<class ValType> void get_unsigned(
+ InputIterator& it, InputIterator end, std::ios_base& iosb,
+ const std::ctype<CharType>& ct,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ switch(peek_char(it, end, ct)) {
+
+ case 'i':
+ get_i(it, end, ct, state, val);
+ break;
+
+ case 'n':
+ get_n(it, end, ct, state, val);
+ break;
+
+ case 'q':
+ case 's':
+ get_q(it, end, ct, state, val);
+ break;
+
+ default:
+ it = std::num_get<CharType, InputIterator>::do_get(
+ it, end, iosb, state, val);
+ if((flags_ & legacy) && val == static_cast<ValType>(1)
+ && peek_char(it, end, ct) == '#')
+ get_one_hash(it, end, ct, state, val);
+ break;
+ }
+ }
+
+ //..........................................................................
+
+ template<class ValType> void get_i(
+ InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ if(!std::numeric_limits<ValType>::has_infinity
+ || (flags_ & trap_infinity)) {
+ state |= std::ios_base::failbit;
+ return;
+ }
+
+ ++it;
+
+ if(!match_string(it, end, ct, "nf")) {
+ state |= std::ios_base::failbit;
+ return;
+ }
+
+ if(peek_char(it, end, ct) != 'i') {
+ val = std::numeric_limits<ValType>::infinity(); // "inf"
+ return;
+ }
+
+ ++it;
+
+ if(!match_string(it, end, ct, "nity")) {
+ state |= std::ios_base::failbit;
+ return;
+ }
+
+ val = std::numeric_limits<ValType>::infinity(); // "infinity"
+ }
+
+ template<class ValType> void get_n(
+ InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ if(!std::numeric_limits<ValType>::has_quiet_NaN
+ || (flags_ & trap_nan)) {
+ state |= std::ios_base::failbit;
+ return;
+ }
+
+ ++it;
+
+ if(!match_string(it, end, ct, "an")) {
+ state |= std::ios_base::failbit;
+ return;
+ }
+
+ switch(peek_char(it, end, ct)) {
+ case 'q':
+ case 's':
+ if(flags_ && legacy)
+ ++it;
+ break; // "nanq", "nans"
+
+ case '(':
+ {
+ ++it;
+ char c;
+ while((c = peek_char(it, end, ct))
+ && c != ')' && c != ' ' && c != '\n' && c != '\t')
+ ++it;
+ if(c != ')') {
+ state |= std::ios_base::failbit;
+ return;
+ }
+ ++it;
+ break; // "nan(...)"
+ }
+
+ default:
+ break; // "nan"
+ }
+
+ val = positive_nan<ValType>();
+ }
+
+ template<class ValType> void get_q(
+ InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ if(!std::numeric_limits<ValType>::has_quiet_NaN
+ || (flags_ & trap_nan) || !(flags_ & legacy)) {
+ state |= std::ios_base::failbit;
+ return;
+ }
+
+ ++it;
+
+ if(!match_string(it, end, ct, "nan")) {
+ state |= std::ios_base::failbit;
+ return;
+ }
+
+ val = positive_nan<ValType>(); // qnan, snan
+ }
+
+ template<class ValType> void get_one_hash(
+ InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ ++it;
+
+ switch(peek_char(it, end, ct)) {
+ case 'i':
+ get_one_hash_i(it, end, ct, state, val);
+ return;
+
+ case 'q':
+ case 's':
+ if(std::numeric_limits<ValType>::has_quiet_NaN
+ && !(flags_ & trap_nan)) {
+ ++it;
+ if(match_string(it, end, ct, "nan")) {
+ // "1.#QNAN", "1.#SNAN"
+ ++it;
+ val = positive_nan<ValType>();
+ return;
+ }
+ }
+ break;
+
+ default:
+ break;
+ }
+
+ state |= std::ios_base::failbit;
+ }
+
+ template<class ValType> void get_one_hash_i(
+ InputIterator& it, InputIterator end, const std::ctype<CharType>& ct,
+ std::ios_base::iostate& state, ValType& val) const
+ {
+ ++it;
+
+ if(peek_char(it, end, ct) == 'n') {
+ ++it;
+ switch(peek_char(it, end, ct)) {
+ case 'f': // "1.#INF"
+ if(std::numeric_limits<ValType>::has_infinity
+ && !(flags_ & trap_infinity)) {
+ ++it;
+ val = std::numeric_limits<ValType>::infinity();
+ return;
+ }
+ break;
+
+ case 'd': // 1.#IND"
+ if(std::numeric_limits<ValType>::has_quiet_NaN
+ && !(flags_ & trap_nan)) {
+ ++it;
+ val = positive_nan<ValType>();
+ return;
+ }
+ break;
+
+ default:
+ break;
+ }
+ }
+
+ state |= std::ios_base::failbit;
+ }
+
+ //..........................................................................
+
+ char peek_char(
+ InputIterator& it, InputIterator end,
+ const std::ctype<CharType>& ct) const
+ {
+ if(it == end) return 0;
+ return ct.narrow(ct.tolower(*it), 0);
+ }
+
+ bool match_string(
+ InputIterator& it, InputIterator end,
+ const std::ctype<CharType>& ct, const char* s) const
+ {
+ while(it != end && *s && *s == ct.narrow(ct.tolower(*it), 0)) {
+ ++s;
+ ++it;
+ }
+ return !*s;
+ }
+
+private:
+ const int flags_;
+};
+
+//------------------------------------------------------------------------------
+
+} // namespace serialization
+} // namespace boost
+
+#ifdef _MSC_VER
+# pragma warning(pop)
+#endif
+
+#endif
=== added file 'extra/floating_point_utilities_v3/boost/math/signbit.hpp'
--- extra/floating_point_utilities_v3/boost/math/signbit.hpp 1970-01-01 00:00:00 +0000
+++ extra/floating_point_utilities_v3/boost/math/signbit.hpp 2014-05-26 08:10:59 +0000
@@ -0,0 +1,86 @@
+// signbit.hpp
+
+#ifndef BOOST_MATH_SIGNBIT_HPP
+#define BOOST_MATH_SIGNBIT_HPP
+
+// Copyright (c) 2006 Johan Rade
+
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include "detail/fp_traits.hpp"
+
+namespace boost {
+namespace math {
+
+//------------------------------------------------------------------------------
+
+template<class T> bool (signbit)(T x)
+{
+ typedef BOOST_DEDUCED_TYPENAME detail::fp_traits<T>::type traits;
+ traits::init();
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= traits::sign;
+ return a != 0;
+}
+
+//------------------------------------------------------------------------------
+
+namespace detail {
+
+ template<class T> T copysign_impl(T x, T y)
+ {
+ typedef BOOST_DEDUCED_TYPENAME fp_traits<T>::type traits;
+ traits::init();
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a &= ~traits::sign;
+
+ BOOST_DEDUCED_TYPENAME traits::bits b;
+ traits::get_bits(y,b);
+ b &= traits::sign;
+
+ traits::set_bits(x,a|b);
+ return x;
+ }
+}
+
+inline float (copysign)(float x, float y) // magnitude of x and sign of y
+{
+ return detail::copysign_impl(x,y);
+}
+
+inline double (copysign)(double x, double y)
+{
+ return detail::copysign_impl(x,y);
+}
+
+inline long double (copysign)(long double x, long double y)
+{
+ return detail::copysign_impl(x,y);
+}
+
+//------------------------------------------------------------------------------
+
+template<class T> T (changesign)(T x)
+{
+ typedef BOOST_DEDUCED_TYPENAME detail::fp_traits<T>::type traits;
+ traits::init();
+
+ BOOST_DEDUCED_TYPENAME traits::bits a;
+ traits::get_bits(x,a);
+ a ^= traits::sign;
+ traits::set_bits(x,a);
+ return x;
+}
+
+//------------------------------------------------------------------------------
+
+} // namespace math
+} // namespace boost
+
+#endif
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2014-04-24 21:37:44 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2014-05-26 23:23:29 +0000
@@ -154,7 +154,7 @@
void measurePressureProfile(double WallUpy, double WallDowny);
double averageSlicePressure(double Y);
double averagePressure();
- double getCell (double X,double Y,double Z);
+ int getCell (double X,double Y,double Z);
double boundaryFlux(unsigned int boundaryId);
vector<Real> averageFluidVelocityOnSphere(unsigned int Id_sph);
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-04-24 21:37:44 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-05-26 23:23:29 +0000
@@ -217,9 +217,9 @@
}
template <class Tesselation>
-double FlowBoundingSphere<Tesselation>::getCell (double X, double Y, double Z)
+int FlowBoundingSphere<Tesselation>::getCell (double X, double Y, double Z)
{
- if (noCache) {cerr<<"Triangulation does not exist. Waht did you do?!"<<endl; return -1;}
+ if (noCache) {cout<<"Triangulation does not exist. Waht did you do?!"<<endl; return -1;}
RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
CellHandle cell = Tri.locate(Point(X,Y,Z));
return cell->info().id;
@@ -348,8 +348,9 @@
}
}
noCache=false;//cache should always be defined after execution of this function
+ }
if (onlyCache) return;
- } else {//use cached values
+// } else {//use cached values
#ifndef parallel_forces
for (FiniteCellsIterator cell = Tri.finite_cells_begin(); cell != cellEnd; cell++) {
for (int yy=0;yy<4;yy++) cell->vertex(yy)->info().forces = cell->vertex(yy)->info().forces + cell->info().unitForceVectors[yy]*cell->info().p();}
@@ -367,7 +368,7 @@
v->info().forces = tf;
}
#endif
- }
+// }
if (debugOut) {
CVector totalForce = nullVect;
for (FiniteVerticesIterator v = Tri.finite_vertices_begin(); v != Tri.finite_vertices_end(); ++v) {
=== modified file 'pkg/dem/CohFrictMat.hpp'
--- pkg/dem/CohFrictMat.hpp 2014-05-19 15:00:11 +0000
+++ pkg/dem/CohFrictMat.hpp 2014-05-26 16:49:23 +0000
@@ -23,7 +23,8 @@
((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 will be elastic."))
+ ((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,0,,"Tensile strength, homogeneous to a pressure"))
((Real,shearCohesion,0,,"Shear strength, homogeneous to a pressure"))
((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."))
=== modified file 'pkg/dem/CohFrictPhys.hpp'
--- pkg/dem/CohFrictPhys.hpp 2014-01-21 15:08:51 +0000
+++ pkg/dem/CohFrictPhys.hpp 2014-05-26 16:49:23 +0000
@@ -23,8 +23,8 @@
((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 to determine the maximum plastic rolling moment."))
- ((Vector3r,maxTwistMoment,Vector3r::Zero(),,"Maximum elastic value for the twisting moment (if zero, plasticity will not be applied). In CohFrictMat a parameter should be added to decide what value should be attributed to this threshold value."))
+ ((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."))
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp 2013-03-06 17:30:45 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp 2014-05-26 16:49:23 +0000
@@ -16,6 +16,9 @@
YADE_PLUGIN((CohesiveFrictionalContactLaw)(Law2_ScGeom6D_CohFrictPhys_CohesionMoment));
CREATE_LOGGER(Law2_ScGeom6D_CohFrictPhys_CohesionMoment);
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::getPlasticDissipation() {return (Real) plasticDissipation;}
+void Law2_ScGeom6D_CohFrictPhys_CohesionMoment::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
+
Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy()
{
Real normEnergy=0;
@@ -41,6 +44,50 @@
return shearEnergy;
}
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::bendingElastEnergy()
+{
+ Real bendingEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+ if (phys) {
+ bendingEnergy += 0.5*(phys->moment_bending.squaredNorm()/phys->kr);
+ }
+ }
+ return bendingEnergy;
+}
+
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::twistElastEnergy()
+{
+ Real twistEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+ if (phys) {
+ twistEnergy += 0.5*(phys->moment_twist.squaredNorm()/phys->ktw);
+ }
+ }
+ return twistEnergy;
+}
+
+Real Law2_ScGeom6D_CohFrictPhys_CohesionMoment::totalElastEnergy()
+{
+ Real totalEnergy=0;
+ FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
+ if(!I->isReal()) continue;
+ CohFrictPhys* phys = YADE_CAST<CohFrictPhys*>(I->phys.get());
+ if (phys) {
+ totalEnergy += 0.5*(phys->normalForce.squaredNorm()/phys->kn);
+ totalEnergy += 0.5*(phys->shearForce.squaredNorm()/phys->ks);
+ totalEnergy += 0.5*(phys->moment_bending.squaredNorm()/phys->kr);
+ totalEnergy += 0.5*(phys->moment_twist.squaredNorm()/phys->ktw);
+ }
+ }
+ return totalEnergy;
+}
+
+
+
void CohesiveFrictionalContactLaw::action()
{
if(!functor) functor=shared_ptr<Law2_ScGeom6D_CohFrictPhys_CohesionMoment>(new Law2_ScGeom6D_CohFrictPhys_CohesionMoment);
@@ -104,17 +151,12 @@
Vector3r trialForce=shearForce;
shearForce *= maxFs;
if (scene->trackEnergy){
- Real dissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
- if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);}
+ Real sheardissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
+ if(sheardissip>0) scene->energy->add(sheardissip,"shearDissip",shearDissipIx,/*reset*/false);}
if (Fn<0) phys->normalForce = Vector3r::Zero();//Vector3r::Zero()
}
//Apply the force
applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
-// Vector3r force = -phys->normalForce-shearForce;
-// scene->forces.addForce(id1,force);
-// scene->forces.addForce(id2,-force);
-// scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
-// scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
/// Moment law ///
if (phys->momentRotationLaw && (!phys->cohesionBroken || always_use_moment_law)) {
@@ -152,22 +194,28 @@
/// Plasticity ///
// limit rolling moment to the plastic value, if required
Real RollMax = phys->maxRollPl*phys->normalForce.norm();
- if (RollMax>0.){ // do we want to apply plasticity?
+ if (RollMax>=0.){ // do we want to apply plasticity?
if (!useIncrementalForm) LOG_WARN("If :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::useIncrementalForm` is false, then plasticity will not be applied correctly (the total formulation would not reproduce irreversibility).");
Real scalarRoll = phys->moment_bending.norm();
if (scalarRoll>RollMax){ // fix maximum rolling moment
Real ratio = RollMax/scalarRoll;
phys->moment_bending *= ratio;
- }
+ if (scene->trackEnergy){
+ Real bendingdissip=((1/phys->kr)*(scalarRoll-RollMax)*RollMax)/*active force*/;
+ if(bendingdissip>0) scene->energy->add(bendingdissip,"bendingDissip",bendingDissipIx,/*reset*/false);}
+ }
}
// limit twisting moment to the plastic value, if required
- Real TwistMax = phys->maxTwistMoment.norm();
- if (TwistMax>0.){ // do we want to apply plasticity?
+ Real TwistMax = phys->maxTwistPl*phys->normalForce.norm();
+ if (TwistMax>=0.){ // do we want to apply plasticity?
if (!useIncrementalForm) LOG_WARN("If :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::useIncrementalForm` is false, then plasticity will not be applied correctly (the total formulation would not reproduce irreversibility).");
Real scalarTwist= phys->moment_twist.norm();
if (scalarTwist>TwistMax){ // fix maximum rolling moment
Real ratio = TwistMax/scalarTwist;
phys->moment_twist *= ratio;
+ if (scene->trackEnergy){
+ Real twistdissip=((1/phys->ktw)*(scalarTwist-TwistMax)*TwistMax)/*active force*/;
+ if(twistdissip>0) scene->energy->add(twistdissip,"twistDissip",twistDissipIx,/*reset*/false);}
}
}
// Apply moments now
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-01-21 15:06:34 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-05-26 16:49:23 +0000
@@ -17,20 +17,32 @@
class Law2_ScGeom6D_CohFrictPhys_CohesionMoment: public LawFunctor{
public:
+ OpenMPAccumulator<Real> plasticDissipation;
Real normElastEnergy();
Real shearElastEnergy();
+ Real bendingElastEnergy();
+ Real twistElastEnergy();
+ Real totalElastEnergy();
+ Real getPlasticDissipation();
+ void initPlasticDissipation(Real initVal=0);
virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_s$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. The maximum value of moments can be defined and takes the form of rolling friction. Cohesive -type moment may also be included in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
((bool,shear_creep,false,,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
((bool,twist_creep,false,,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
+ ((bool,traceEnergy,false,,"Define the total energy dissipated in plastic slips at all contacts. This will trace only plastic energy in this law, see O.trackEnergy for a more complete energies tracing"))
((bool,useIncrementalForm,false,,"use the incremental formulation to compute bending and twisting moments. Creep on the twisting moment is not included in such a case."))
- ((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
+ ((int,shearDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for shear dissipation (with O.trackEnergy)"))
+ ((int,bendingDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for bending dissipation (with O.trackEnergy)"))
+ ((int,twistDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for twist dissipation (with O.trackEnergy)"))
((Real,creep_viscosity,1,,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_CohFrictMat_CohFrictMat_CohFrictPhys."))
,,
.def("normElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
.def("shearElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
+ .def("bendingElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::bendingElastEnergy,"Compute bending elastic energy.")
+ .def("twistElastEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::twistElastEnergy,"Compute twist elastic energy.")
+ .def("elasticEnergy",&Law2_ScGeom6D_CohFrictPhys_CohesionMoment::totalElastEnergy,"Compute total elastic energy.")
);
FUNCTOR2D(ScGeom6D,CohFrictPhys);
DECLARE_LOGGER;
=== modified file 'pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp 2014-01-24 17:15:20 +0000
+++ pkg/dem/Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.cpp 2014-05-26 16:49:23 +0000
@@ -51,9 +51,6 @@
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;
- // Jean-Patrick Plassiard, Noura Belhaine, Frederic
- // Victor Donze, "Calibration procedure for spherical
- // discrete elements using a local moemnt law".
contactPhysics->kr = Da*Db*Ks*AlphaKr;
contactPhysics->ktw = Da*Db*Ks*AlphaKtw;
contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
@@ -69,8 +66,8 @@
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);
- //contactPhysics->elasticRollingLimit = elasticRollingLimit;
}
else {// !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
CohFrictPhys* contactPhysics = YADE_CAST<CohFrictPhys*>(interaction->phys.get());
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.cpp'
--- pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-05-23 13:05:19 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-05-28 08:47:32 +0000
@@ -41,7 +41,7 @@
}
if ( smoothJoint && phys->isOnJoint ) {
- if ( phys->more || ( phys->jointCumulativeSliding > (2*min(geom->radius1,geom->radius2)) ) ) {
+ if ( phys->more || ( phys->jointCumulativeSliding > (2*min(geom->radius1,geom->radius2)) ) ) {
scene->interactions->requestErase(contact); return;
} else {
D = phys->initD - abs((b1->state->pos - b2->state->pos).dot(phys->jointNormal));
@@ -51,8 +51,11 @@
}
/* Determination of interaction */
- if (D < 0) { //spheres do not touch
- if ( !phys->isCohesive ) { scene->interactions->requestErase(contact); return; }
+ if (D < 0) { //spheres do not "touch" !
+ if ( !phys->isCohesive )
+ {
+ scene->interactions->requestErase(contact); return;
+ }
if ( phys->isCohesive && (phys->FnMax>0) && (abs(D)>Dtensile) ) {
@@ -138,9 +141,12 @@
cracksFileExist=true;
// delete contact if in tension, set the contact properties to friction if in compression
- if ( D < 0 ) {
+ if ( D < 0 )
+ {
scene->interactions->requestErase(contact); return;
- } else {
+ }
+ else
+ {
phys->FnMax = 0;
phys->FsMax = 0;
phys->isCohesive=false;
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.hpp'
--- pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-03-14 08:46:30 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-05-28 08:47:32 +0000
@@ -98,7 +98,7 @@
virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
FUNCTOR2D(ScGeom,JCFpmPhys);
- YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM,LawFunctor,"Interaction law for cohesive frictional material, e.g. rock, possibly presenting joint surfaces. Joint surfaces can be defined in a preprocessing phase through .stl meshes (see ref for details of the procedure), and can be mechanically described with a smooth contact logic [Ivars2011]_ (implemented in Yade in [Scholtes2012]_).",
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM,LawFunctor,"Interaction law for cohesive frictional material, e.g. rock, possibly presenting joint surfaces, that can be mechanically described with a smooth contact logic [Ivars2011]_ (implemented in Yade in [Scholtes2012]_). See examples/jointedCohesiveFrictionalPM for script examples. Joint surface definitions (through stl meshes or direct definition with gts module) are illustrated there.",
((string,Key,"",,"string specifying the name of saved file 'cracks___.txt', when :yref:`recordCracks<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>` is true."))
((bool,cracksFileExist,false,,"if true (and if :yref:`recordCracks<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.recordCracks>`), data are appended to an existing 'cracksKey' text file; otherwise its content is reset."))
((bool,smoothJoint,false,,"if true, interactions of particles belonging to joint surface (:yref:`JCFpmPhys.isOnJoint`) are handled according to a smooth contact logic [Ivars2011]_, [Scholtes2012]_."))
=== modified file 'pkg/dem/KinemSimpleShearBox.cpp'
--- pkg/dem/KinemSimpleShearBox.cpp 2013-05-29 09:48:51 +0000
+++ pkg/dem/KinemSimpleShearBox.cpp 2014-05-28 08:55:24 +0000
@@ -55,24 +55,16 @@
Real Ysup = topbox->state->pos.y();
Real Ylat = leftbox->state->pos.y();
-// Changes in vertical and horizontal position :
-
-
- topbox->state->pos += Vector3r(dX,dY,0);
-
- leftbox->state->pos += Vector3r(dX/2.0,dY/2.0,0);
- rightbox->state->pos += Vector3r(dX/2.0,dY/2.0,0);
- if(LOG) cout << "dY reellemt applique :" << dY << endl;
- if(LOG) cout << "qui nous a emmene en : y = " <<(topbox->state->pos).y() << endl;
-
- Real Ysup_mod = topbox->state->pos.y();
- Real Ylat_mod = leftbox->state->pos.y();
-
-// with the corresponding velocities :
+// Changes in vertical and horizontal velocities :
topbox->state->vel = Vector3r(dX/dt,dY/dt,0);
leftbox->state->vel = Vector3r(dX/(2.0 * dt),dY/(2.0 * dt),0);
rightbox->state->vel = Vector3r(dX/(2.0 * dt),dY/(2.0*dt),0);
+ if(LOG) cout << "dY that will be applied by NewtonIntegrator :" << dY << endl;
+
+ Real Ysup_mod = Ysup + dY;
+ Real Ylat_mod = Ylat + dY;
+
computeAlpha();
// Then computation of the angle of the rotation,dalpha, to be done :
if (alpha == Mathr::PI/2.0) // Case of the very beginning
@@ -87,11 +79,8 @@
Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));
-// Rotation is applied : orientation and angular velocity of plates are modified.
- leftbox->state->ori = qcorr*leftbox->state->ori;
+// Rotation is applied through velocities (and NewtonIntegrator)
leftbox->state->angVel = Vector3r(0,0,1)*dalpha/dt;
-
- rightbox->state->ori = qcorr*rightbox->state->ori;
rightbox->state->angVel = Vector3r(0,0,1)*dalpha/dt;
}
=== modified file 'pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp'
--- pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 2013-06-25 08:02:13 +0000
+++ pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 2014-05-26 09:09:48 +0000
@@ -15,8 +15,7 @@
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,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. The maximum value of moments can be defined and takes the form of rolling friction. Cohesive -type moment may also be included in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
- ((bool,useIncrementalForm,false,,"use the incremental formulation to compute bending and twisting moments. Creep on the twisting moment is not included in such a case."))
+ 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.")
=== modified file 'pkg/dem/NormalInelasticityLaw.hpp'
--- pkg/dem/NormalInelasticityLaw.hpp 2011-05-11 09:40:28 +0000
+++ pkg/dem/NormalInelasticityLaw.hpp 2014-05-28 08:17:25 +0000
@@ -33,7 +33,7 @@
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/normalInelasticityTest.py",
+ "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"))
,
=== modified file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 2014-05-08 05:57:04 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-05-27 11:24:16 +0000
@@ -118,8 +118,12 @@
NLiqBridg -= 1;
}
#ifdef YADE_LIQMIGRATION
- const intReal B1={id1, phys.Vb/2.0};
- const intReal B2={id2, phys.Vb/2.0};
+ if (phys.Vb > 0.0 and ((phys.Vf1+phys.Vf2) == 0.0)) {
+ phys.Vf1 = phys.Vb/2.0;
+ phys.Vf2 = phys.Vb/2.0;
+ }
+ const intReal B1={id1, phys.Vf1};
+ const intReal B2={id2, phys.Vf2};
scene->delIntrs.push_back(B1);
scene->delIntrs.push_back(B2);
#endif
@@ -411,6 +415,13 @@
addBodyMapReal(bodyUpdateLiquid, id2, -Vf2);
Vb->Vb = Vrup;
+ if (particleconserve) {
+ Vb->Vf1 = Vf1;
+ Vb->Vf2 = Vf2;
+ } else {
+ Vb->Vf1 = Vrup/2.0;
+ Vb->Vf2 = Vrup/2.0;
+ }
}
scene->addIntrs.clear();
@@ -462,8 +473,19 @@
if(!((*it).second) or !(((*it).second)->isReal())) continue;
ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
if (physT->Vb<physT->Vmax) {
- liqVolShr += (physT->Vmax - physT->Vb)*FillLevel;
- physT->Vb += (physT->Vmax - physT->Vb)*FillLevel;
+ const Real addVolLiq = (physT->Vmax - physT->Vb)*FillLevel;
+ liqVolShr += addVolLiq;
+ physT->Vb += addVolLiq;
+ if (particleconserve) {
+ if (((*it).second)->getId1() == (*it).first) {
+ physT->Vf1+=addVolLiq;
+ } else if (((*it).second)->getId2() == (*it).first) {
+ physT->Vf2+=addVolLiq;
+ }
+ } else {
+ physT->Vf1+=addVolLiq/2.0;
+ physT->Vf2+=addVolLiq/2.0;
+ }
}
}
return;
=== modified file 'pkg/dem/ViscoelasticCapillarPM.hpp'
--- pkg/dem/ViscoelasticCapillarPM.hpp 2014-05-23 13:20:43 +0000
+++ pkg/dem/ViscoelasticCapillarPM.hpp 2014-05-27 14:20:39 +0000
@@ -1,5 +1,6 @@
#include"ViscoelasticPM.hpp"
#include <boost/unordered_map.hpp>
+#include<yade/core/PartialEngine.hpp>
class ViscElCapMat : public ViscElMat {
public:
@@ -35,6 +36,8 @@
((CapType,CapillarType,None_Capillar,,"Different types of capillar interaction: Willett_numeric, Willett_analytic, Weigert, Rabinovich, Lambert, Soulie"))
#ifdef YADE_LIQMIGRATION
((Real,Vmax,0.0,,"Maximal liquid bridge volume [m^3]"))
+ ((Real,Vf1,0.0,, "Liquid which will be returned to the 1st body after rupture [m^3]"))
+ ((Real,Vf2,0.0,, "Liquid which will be returned to the 2nd body after rupture [m^3]"))
#endif
,
createIndex();
@@ -95,6 +98,7 @@
((Real,liqVolRup,0.,, "Liquid volume (integral value), which has been freed after rupture occured, [m^3]."))
((Real,liqVolShr,0.,, "Liquid volume (integral value), which has been shared among of contacts, [m^3]."))
((Real,vMaxCoef,0.03,, "Coefficient for vMax, [-]."))
+ ((bool,particleconserve,false,, "If True, the particle will have the same liquid volume during simulation e.g. liquid will not migrate [false]."))
,/* ctor */
,/* py */
.def("totalLiq",&LiqControl::totalLiqVol,(boost::python::arg("mask")=0),"Return total volume of water in simulation.")
=== modified file 'pkg/pfv/DFNFlow.cpp'
--- pkg/pfv/DFNFlow.cpp 2014-04-03 16:44:33 +0000
+++ pkg/pfv/DFNFlow.cpp 2014-05-27 10:33:52 +0000
@@ -38,9 +38,11 @@
void trickPermeability();
void trickPermeability (RTriangulation::Facet_circulator& facet,Real somethingBig);
void trickPermeability (RTriangulation::Finite_edges_iterator& edge,Real somethingBig);
+ void setPositionsBuffer(bool current);
YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DFNFlowEngine,DFNFlowEngineT,"documentation here",
((Real, myNewAttribute, 0,,"useless example"))
+ ((bool, updatePositions, false,,"update particles positions when rebuilding the mesh (experimental)"))
,/*DFNFlowEngineT()*/,
,
)
@@ -48,6 +50,26 @@
};
REGISTER_SERIALIZABLE(DFNFlowEngine);
YADE_PLUGIN((DFNFlowEngine));
+//In this version, we never update positions when !updatePositions, i.e. keep triangulating the same positions
+void DFNFlowEngine::setPositionsBuffer(bool current)
+{
+ vector<posData>& buffer = current? positionBufferCurrent : positionBufferParallel;
+ if (!updatePositions && buffer.size()>0) return;
+ buffer.clear();
+ buffer.resize(scene->bodies->size());
+ shared_ptr<Sphere> sph ( new Sphere );
+ const int Sph_Index = sph->getClassIndexStatic();
+ FOREACH ( const shared_ptr<Body>& b, *scene->bodies ) {
+ if (!b || ignoredBody==b->getId()) continue;
+ posData& dat = buffer[b->getId()];
+ dat.id=b->getId();
+ dat.pos=b->state->pos;
+ dat.isSphere= (b->shape->getClassIndex() == Sph_Index);
+ if (dat.isSphere) dat.radius = YADE_CAST<Sphere*>(b->shape.get())->radius;
+ dat.exists=true;
+ }
+}
+
void DFNFlowEngine::trickPermeability (RTriangulation::Facet_circulator& facet, Real somethingBig)
{
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-05-23 13:05:19 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-05-27 07:50:01 +0000
@@ -79,7 +79,7 @@
typedef typename RTriangulation::Triangulation_data_structure::Cell::Info CellInfo;
typedef typename RTriangulation::Triangulation_data_structure::Vertex::Info VertexInfo;
- protected:
+// protected:
shared_ptr<FlowSolver> solver;
shared_ptr<FlowSolver> backgroundSolver;
volatile bool backgroundCompleted;
@@ -88,7 +88,7 @@
vector<posData> positionBufferCurrent;//reflect last known positions before we start computations
vector<posData> positionBufferParallel;//keep the positions from a given step for multithread factorization
//copy positions in a buffer for faster and/or parallel access
- void setPositionsBuffer(bool current);
+ virtual void setPositionsBuffer(bool current);
virtual void trickPermeability() {};
public :
=== modified file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 2014-05-07 10:28:20 +0000
+++ pkg/pfv/FlowEngine.ipp 2014-05-26 23:23:29 +0000
@@ -249,10 +249,7 @@
void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( double pZero, Solver& flow )
{
if (first) flow.currentTes=0;
- else {
- flow.currentTes=!flow.currentTes;
- if (debug) cout << "--------RETRIANGULATION-----------" << endl;
- }
+ else { flow.currentTes=!flow.currentTes; if (debug) cout << "--------RETRIANGULATION-----------" << endl;}
flow.resetNetwork();
initSolver(flow);
=== modified file 'py/config.py.in'
--- py/config.py.in 2014-02-28 08:22:32 +0000
+++ py/config.py.in 2014-05-27 11:24:16 +0000
@@ -14,7 +14,10 @@
libDir=os.path.abspath(prefix+'/'+libPATH+'/yade${SUFFIX}')
confDir=os.environ['HOME']+'/.yade${SUFFIX}'
libstdcxx='${libstdcxx}'
-features='${features}'.split(' ')
+
+features='${CONFIGURED_FEATS}'.split(' ')
+if (features[0]==''): features=features[1:]
+
revision='${realVersion}'
version='${version}'
sourceRoot='${sourceRoot}'
=== modified file 'py/ymport.py'
--- py/ymport.py 2014-05-20 06:35:16 +0000
+++ py/ymport.py 2014-05-26 07:53:31 +0000
@@ -84,6 +84,10 @@
if (len(curClump)<>0):
ret.append(O.bodies.appendClumped(curClump,discretization=discretization))
+
+ # Set the mask to a clump the same as the first member of it
+ for i in range(len(ret)):
+ O.bodies[ret[i][0]].mask = O.bodies[ret[i][1][0]].mask
return ret
def text(fileName,shift=Vector3.Zero,scale=1.0,**kw):
=== modified file 'scripts/checks-and-tests/checks/DEM-PFV-check.py'
--- scripts/checks-and-tests/checks/DEM-PFV-check.py 2014-01-23 18:10:35 +0000
+++ scripts/checks-and-tests/checks/DEM-PFV-check.py 2014-05-27 11:29:38 +0000
@@ -2,11 +2,8 @@
# Here, we are testing bulk modulus, then permeability, then the consolidation of a specimen.
# the test is based on examples/FluidCouplingPFV/oedometer.py, only slightly simplified and using less particles
-try: FlowEngine
-except NameError:
- print "skip DEM-PFV check, FlowEngine not available"
-else:
+if ('PFVflow' in features):
errors=0
tolerance=0.01
@@ -143,4 +140,6 @@
errors+=1
if (errors):
- resultStatus +=1 #Test is failed
\ No newline at end of file
+ resultStatus +=1 #Test is failed
+else:
+ print "skip DEM-PFV check, FlowEngine not available"
=== added file 'scripts/checks-and-tests/checks/checkLiquidMigration.py'
--- scripts/checks-and-tests/checks/checkLiquidMigration.py 1970-01-01 00:00:00 +0000
+++ scripts/checks-and-tests/checks/checkLiquidMigration.py 2014-05-27 11:24:16 +0000
@@ -0,0 +1,104 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+# The model checks liquid migration model if it is enabled during compilation
+from yade import utils, plot
+
+if ('LIQMIGRATION' in features):
+ o = Omega()
+ fr = 0.5;rho=2000
+ tc = 0.001; en = 0.7; et = 0.7;
+ o.dt = 1.0
+
+
+ r1 = 1.0
+ r2 = 1.0
+ Gamma = 20.6*1e-3
+ Theta = 0
+ VB = 74.2*1e-12
+
+ tolerance = 1e-6
+
+
+ CapillarType = "Lambert"
+
+ mat1 = O.materials.append(ViscElCapMat(frictionAngle=fr,density=rho,Vb=VB,gamma=Gamma,theta=Theta,Capillar=True,CapillarType=CapillarType,tc=tc,en=en,et=et))
+
+ d = 1.1
+ id0 = O.bodies.append(sphere(center=[0,0,0],radius=r1,material=mat1,fixed=True, color=[1,0,0]))
+
+ id1 = O.bodies.append(sphere(center=[0,0,(r1+r2)*d],radius=r2,material=mat1,fixed=True, color=[0,1,0]))
+ id2 = O.bodies.append(sphere(center=[0,0,-(r1+r2)*d],radius=r2,material=mat1,fixed=True, color=[0,1,0]))
+
+
+ O.bodies[id0].Vf = 0.3e-1
+ O.bodies[id0].Vmin = 0.1e-1
+
+ O.bodies[id1].Vf = 0.4e-1
+ O.bodies[id1].Vmin = 0.1e-1
+
+ O.bodies[id2].Vf = 0.5e-1
+ O.bodies[id2].Vmin = 0.1e-1
+
+ vel = -0.15
+ O.bodies[id1].state.vel=[0,0,vel]
+ O.bodies[id2].state.vel=[0,0,-vel]
+
+ o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=(r1+r2)*5.0),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom()],
+ [Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],
+ [Law2_ScGeom_ViscElCapPhys_Basic()],
+ ),
+ LiqControl(particleconserve=True,label='liqcontrol'),
+ NewtonIntegrator(damping=0,gravity=[0,0,0]),
+ PyRunner(command='showData()',iterPeriod=1,dead=True),
+ ]
+
+ def showData():
+ print "Step %d"%O.iter
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id0, O.bodies[id0].Vf, O.bodies[id0].Vmin)
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id1, O.bodies[id1].Vf, O.bodies[id1].Vmin)
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id2, O.bodies[id2].Vf, O.bodies[id2].Vmin)
+ try:
+ print "Interaction[%d, %d].Vb=%s"%(id0, id1, O.interactions[id0,id1].phys.Vb)
+ except:
+ pass
+
+ try:
+ print "Interaction[%d, %d].Vb=%s"%(id0, id2, O.interactions[id0,id2].phys.Vb)
+ except:
+ pass
+ print
+
+ def switchVel():
+ O.bodies[id1].state.vel=-O.bodies[id1].state.vel
+ O.bodies[id2].state.vel=-O.bodies[id2].state.vel
+
+ resultStatus = 0
+ O.run(3, True)
+ if ((abs((O.interactions[id0,id1].phys.Vb - 0.03)/0.03) > tolerance) or
+ (abs((O.interactions[id0,id1].phys.Vb - 0.03)/0.03) > tolerance)):
+ resultStatus += 1
+
+ switchVel()
+ O.run(5, True)
+ if ((abs((O.bodies[id0].Vf - 0.03)/0.03) > tolerance) or
+ (abs((O.bodies[id1].Vf - 0.04)/0.04) > tolerance) or
+ (abs((O.bodies[id2].Vf - 0.05)/0.05) > tolerance)):
+ resultStatus += 1
+
+ liqcontrol.particleconserve=False
+ switchVel()
+ O.run(5, True)
+ switchVel()
+ O.run(5, True)
+ if ((abs((O.bodies[id0].Vf - 0.0465)/0.0465) > tolerance) or
+ (abs((O.bodies[id1].Vf - 0.0325)/0.0325) > tolerance) or
+ (abs((O.bodies[id2].Vf - 0.041)/0.041) > tolerance)):
+ resultStatus += 1
+else:
+ print "This checkLiquidMigration.py cannot be executed because LIQMIGRATION is disables"
+
