yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11436
[Branch ~yade-pkg/yade/git-trunk] Rev 3432: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Christian Jakob (jakob-ifgt)
Jan Stránský (honzik)
Jérôme Duriez (jduriez)
Klaus Thoeni (klaus.thoeni)
------------------------------------------------------------
revno: 3432 [merge]
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Mon 2014-06-23 17:30:42 +0200
message:
Merge github.com:yade/trunk into chaoUnsat
added:
doc/sphinx/ipython_directive200.py
pkg/pfv/FlowEngine.hpp.in
pkg/pfv/FlowEngine.ipp.in
scripts/update_names.sh
modified:
CMakeLists.txt
core/Body.cpp
core/Body.hpp
core/State.hpp
doc/sphinx/conf.py
doc/sphinx/installation.rst
examples/capillary/liquidmigration/showcase.py
examples/concrete/uniax.py
examples/gts-horse/gts-random-pack-obb.py
examples/test/performance/checkPerf.py
gui/qt4/SerializableEditor.py
lib/base/Math.cpp
lib/base/Math.hpp
lib/triangulation/FlowBoundingSphere.ipp
pkg/common/InteractionLoop.cpp
pkg/dem/CapillaryTriaxialTest.cpp
pkg/dem/CohesiveTriaxialTest.cpp
pkg/dem/FrictViscoPM.hpp
pkg/dem/JointedCohesiveFrictionalPM.cpp
pkg/dem/JointedCohesiveFrictionalPM.hpp
pkg/dem/SimpleShear.cpp
pkg/dem/SpheresFactory.cpp
pkg/dem/TesselationWrapper.cpp
pkg/dem/TriaxialTest.cpp
pkg/dem/VTKRecorder.cpp
pkg/dem/VTKRecorder.hpp
pkg/dem/ViscoelasticCapillarPM.cpp
pkg/dem/ViscoelasticCapillarPM.hpp
pkg/pfv/DFNFlow.cpp
pkg/pfv/DummyFlowEngine.cpp
pkg/pfv/FlowEngine.cpp
pkg/pfv/FlowEngine.hpp
pkg/pfv/FlowEngine.ipp
pkg/pfv/PeriodicFlowEngine.cpp
pkg/pfv/SoluteFlowEngine.cpp
py/bodiesHandling.py
py/geom.py
py/wrapper/customConverters.cpp
py/ymport.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-30 11:58:01 +0000
+++ CMakeLists.txt 2014-06-20 15:28:21 +0000
@@ -17,9 +17,10 @@
# 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)
# CHUNKSIZE: set >1, if you want several sources to be compiled at once. Increases compilation speed and RAM-consumption during it (1 by default).
-
+# VECTORIZE: enables vectorization and alignment in Eigen3 library, experimental (OFF by default)
project(Yade C CXX)
cmake_minimum_required(VERSION 2.8)
@@ -125,13 +126,27 @@
OPTION(ENABLE_SPH "Enable SPH" OFF)
OPTION(ENABLE_LBMFLOW "Enable LBM engine (very experimental)" OFF)
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)
#===========================================================
# Use Eigen3 by default
IF (EIGEN3_FOUND)
INCLUDE_DIRECTORIES(${EIGEN3_INCLUDE_DIR})
- MESSAGE(STATUS "Found Eigen3")
+ MESSAGE(STATUS "Found Eigen3, version: ${EIGEN3_VERSION}")
+
+ # Minimal required version 3.2.1
+ IF ((${EIGEN3_MAJOR_VERSION} LESS 2) OR ((${EIGEN3_MAJOR_VERSION} EQUAL 2) AND (${EIGEN3_MINOR_VERSION} LESS 1)))
+ 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")
+ ELSE (NOT VECTORIZE)
+ MESSAGE(STATUS "Enable vectorization")
+ ENDIF (NOT VECTORIZE)
+
ENDIF(EIGEN3_FOUND)
#===========================================================
INCLUDE_DIRECTORIES(${BZIP2_INCLUDE_DIR})
@@ -343,6 +358,18 @@
ELSE(ENABLE_LBMFLOW)
SET(DISABLED_FEATS "${DISABLED_FEATS} LBMFLOW")
ENDIF(ENABLE_LBMFLOW)
+#===============================================
+IF(ENABLE_MASK_ARBITRARY)
+ ADD_DEFINITIONS("-DYADE_MASK_ARBITRARY")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} MASK_ARBITRARY")
+ IF(NOT MASK_ARBITRARY_SIZE)
+ SET(MASK_ARBITRARY_SIZE "256")
+ ENDIF(NOT MASK_ARBITRARY_SIZE)
+ ADD_DEFINITIONS(-DYADE_MASK_ARBITRARY_SIZE=${MASK_ARBITRARY_SIZE})
+ MESSAGE("MASK_ARBITRARY_SIZE = ${MASK_ARBITRARY_SIZE}")
+ELSE(ENABLE_MASK_ARBITRARY)
+ SET(DISABLED_FEATS "${DISABLED_FEATS} MASK_ARBITRARY")
+ENDIF(ENABLE_MASK_ARBITRARY)
#===========================================================
@@ -472,6 +499,17 @@
CONFIGURE_FILE(py/__init__.py.in "${CMAKE_BINARY_DIR}/__init__.py")
#===========================================================
+# Create header files for PFV from FlowEngine.hpp.in-template.
+# All @TEMPLATE_FLOW_NAME@ are replacing by a given names
+
+SET (TEMPLATE_FLOW_NAMES DFNFlowEngineT DummyFlowEngineT FlowEngineT FlowEngine_PeriodicInfo SoluteFlowEngineT)
+FOREACH(TF ${TEMPLATE_FLOW_NAMES})
+ SET (TEMPLATE_FLOW_NAME ${TF})
+ CONFIGURE_FILE(pkg/pfv/FlowEngine.hpp.in "${CMAKE_BINARY_DIR}/pkg/pfv/FlowEngine_${TF}.hpp" @ONLY)
+ CONFIGURE_FILE(pkg/pfv/FlowEngine.ipp.in "${CMAKE_BINARY_DIR}/pkg/pfv/FlowEngine_${TF}.ipp" @ONLY)
+ENDFOREACH(TF)
+INCLUDE_DIRECTORIES("${CMAKE_BINARY_DIR}/pkg/pfv/")
+#===========================================================
INSTALL(PROGRAMS "${CMAKE_BINARY_DIR}/bins/yade${SUFFIX}-batch" DESTINATION ${YADE_EXEC_PATH}/)
INSTALL(PROGRAMS "${CMAKE_BINARY_DIR}/bins/yade${SUFFIX}" DESTINATION ${YADE_EXEC_PATH}/)
=== modified file 'core/Body.cpp'
--- core/Body.cpp 2014-05-23 13:05:19 +0000
+++ core/Body.cpp 2014-06-17 10:18:00 +0000
@@ -31,3 +31,13 @@
return intrSize;
}
+
+bool Body::maskOk(int mask) const { return (mask==0 || ((groupMask & mask) != 0)); }
+bool Body::maskCompatible(int mask) const { return (groupMask & mask) != 0; }
+#ifdef YADE_MASK_ARBITRARY
+bool Body::maskOk(const mask_t& mask) const { return (mask==0 || ((groupMask & mask) != 0)); }
+bool Body::maskCompatible(const mask_t& mask) const { return (groupMask & mask) != 0; }
+#endif
+
+
+
=== modified file 'core/Body.hpp'
--- core/Body.hpp 2014-06-05 13:19:44 +0000
+++ core/Body.hpp 2014-06-17 10:18:00 +0000
@@ -21,6 +21,8 @@
#include<yade/lib/multimethods/Indexable.hpp>
+
+
class Scene;
class Interaction;
@@ -31,12 +33,6 @@
// internal structure to hold some interaction of a body; used by InteractionContainer;
typedef std::map<Body::id_t, shared_ptr<Interaction> > MapId2IntrT;
// groupMask type
-#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
- typedef boost::python::long_ groupMask_t;
-#else
- typedef int groupMask_t;
-#endif
-
// bits for Body::flags
enum { FLAG_BOUNDED=1, FLAG_ASPHERICAL=2 }; /* add powers of 2 as needed */
@@ -75,12 +71,12 @@
Body::id_t getId() const {return id;};
unsigned int coordNumber(); // Number of neighboring particles
- Body::groupMask_t getGroupMask() const {return groupMask; };
- bool maskOk(int mask) const {return (mask==0 || (bool)(groupMask & mask));}
- bool maskCompatible(int mask) const { return (bool)(groupMask & mask); }
-#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
- bool maskOk(const boost::python::long_& mask) const {return (mask==0 || (bool)(groupMask & mask));}
- bool maskCompatible(const boost::python::long_& mask) const { return (bool)(groupMask & mask); }
+ mask_t getGroupMask() const {return groupMask; };
+ bool maskOk(int mask) const;
+ bool maskCompatible(int mask) const;
+#ifdef YADE_MASK_ARBITRARY
+ bool maskOk(const mask_t& mask) const;
+ bool maskCompatible(const mask_t& mask) const;
#endif
// only BodyContainer can set the id of a body
@@ -89,7 +85,7 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Body,Serializable,"A particle, basic element of simulation; interacts with other bodies.",
((Body::id_t,id,Body::ID_NONE,Attr::readonly,"Unique id of this body."))
- ((Body::groupMask_t,groupMask,1,,"Bitmask for determining interactions."))
+ ((mask_t,groupMask,1,,"Bitmask for determining interactions."))
((int,flags,FLAG_BOUNDED,Attr::readonly,"Bits of various body-related flags. *Do not access directly*. In c++, use isDynamic/setDynamic, isBounded/setBounded, isAspherical/setAspherical. In python, use :yref:`Body.dynamic`, :yref:`Body.bounded`, :yref:`Body.aspherical`."))
((shared_ptr<Material>,material,,,":yref:`Material` instance associated with this body."))
@@ -119,7 +115,7 @@
.add_property("dynamic",&Body::isDynamic,&Body::setDynamic,"Whether this body will be moved by forces. (In c++, use ``Body::isDynamic``/``Body::setDynamic``) :ydefault:`true`")
.add_property("bounded",&Body::isBounded,&Body::setBounded,"Whether this body should have :yref:`Body.bound` created. Note that bodies without a :yref:`bound <Body.bound>` do not participate in collision detection. (In c++, use ``Body::isBounded``/``Body::setBounded``) :ydefault:`true`")
.add_property("aspherical",&Body::isAspherical,&Body::setAspherical,"Whether this body has different inertia along principal axes; :yref:`NewtonIntegrator` makes use of this flag to call rotation integration routine for aspherical bodies, which is more expensive. :ydefault:`false`")
- .def_readwrite("mask",&Body::groupMask,"Shorthand for :yref:`Body::groupMask`")
+ .add_property("mask",boost::python::make_getter(&Body::groupMask,boost::python::return_value_policy<boost::python::return_by_value>()),boost::python::make_setter(&Body::groupMask,boost::python::return_value_policy<boost::python::return_by_value>()),"Shorthand for :yref:`Body::groupMask`")
.add_property("isStandalone",&Body::isStandalone,"True if this body is neither clump, nor clump member; false otherwise.")
.add_property("isClumpMember",&Body::isClumpMember,"True if this body is clump member, false otherwise.")
.add_property("isClump",&Body::isClump,"True if this body is clump itself, false otherwise.")
@@ -135,16 +131,3 @@
);
};
REGISTER_SERIALIZABLE(Body);
-
-
-#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
-namespace boost { namespace serialization {
-template<class Archive>
-void serialize(Archive & ar, boost::python::long_ & l, const unsigned int version){
- namespace bp = boost::python;
- std::string value = bp::extract<std::string>(bp::str(l));
- ar & BOOST_SERIALIZATION_NVP(value);
- l = bp::long_(value);
-}
-}} // namespace boost::serialization
-#endif
=== modified file 'core/State.hpp'
--- core/State.hpp 2014-03-27 14:52:07 +0000
+++ core/State.hpp 2014-06-12 16:54:31 +0000
@@ -46,9 +46,9 @@
// python access functions: pos and ori are references to inside Se3r and cannot be pointed to directly
Vector3r pos_get() const {return pos;}
- void pos_set(const Vector3r& p) {pos=p;}
+ void pos_set(const Vector3r p) {pos=p;}
Quaternionr ori_get() const {return ori; }
- void ori_set(const Quaternionr& o){ori=o;}
+ void ori_set(const Quaternionr o){ori=o;}
YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(State,Serializable,"State of a body (spatial configuration, internal variables).",
((Se3r,se3,Se3r(Vector3r::Zero(),Quaternionr::Identity()),,"Position and orientation as one object."))
=== modified file 'doc/sphinx/conf.py'
--- doc/sphinx/conf.py 2013-08-14 08:01:55 +0000
+++ doc/sphinx/conf.py 2014-06-07 08:15:18 +0000
@@ -364,8 +364,10 @@
else:
if 12<=yade.runtime.ipython_version<13:
import ipython_directive012 as id
- else:
+ elif 13<=yade.runtime.ipython_version<200:
import ipython_directive013 as id
+ else:
+ import ipython_directive200 as id
#The next four lines are for compatibility with IPython 0.13.1
ipython_rgxin =re.compile(r'(?:In |Yade )\[(\d+)\]:\s?(.*)\s*')
@@ -412,8 +414,10 @@
else:
if 12<=yade.runtime.ipython_version<13:
extensions.append('ipython_directive012')
- else:
+ elif 13<=yade.runtime.ipython_version<200:
extensions.append('ipython_directive013')
+ else:
+ extensions.append('ipython_directive200')
# the sidebar extension
if False:
=== modified file 'doc/sphinx/installation.rst'
--- doc/sphinx/installation.rst 2014-05-14 12:31:59 +0000
+++ doc/sphinx/installation.rst 2014-06-19 15:35:52 +0000
@@ -115,7 +115,7 @@
* `libQGLViewer <http://www.libqglviewer.com>`_
* `python <http://www.python.org>`_, `numpy <http://numpy.scipy.org>`_, `ipython <http://ipython.scipy.org>`_
* `matplotlib <http://matplotlib.sf.net>`_
-* `eigen3 <http://eigen.tuxfamily.org>`_ algebra library
+* `eigen3 <http://eigen.tuxfamily.org>`_ algebra library (minimal required version 3.2.1)
* `gdb <http://www.gnu.org/software/gdb>`_ debugger
* `sqlite3 <http://www.sqlite.org>`_ database engine
* `Loki <http://loki-lib.sf.net>`_ library
=== added file 'doc/sphinx/ipython_directive200.py'
--- doc/sphinx/ipython_directive200.py 1970-01-01 00:00:00 +0000
+++ doc/sphinx/ipython_directive200.py 2014-06-07 08:15:18 +0000
@@ -0,0 +1,1154 @@
+# -*- coding: utf-8 -*-
+"""
+Sphinx directive to support embedded IPython code.
+
+This directive allows pasting of entire interactive IPython sessions, prompts
+and all, and their code will actually get re-executed at doc build time, with
+all prompts renumbered sequentially. It also allows you to input code as a pure
+python input by giving the argument python to the directive. The output looks
+like an interactive ipython section.
+
+To enable this directive, simply list it in your Sphinx ``conf.py`` file
+(making sure the directory where you placed it is visible to sphinx, as is
+needed for all Sphinx directives). For example, to enable syntax highlighting
+and the IPython directive::
+
+ extensions = ['IPython.sphinxext.ipython_console_highlighting',
+ 'IPython.sphinxext.ipython_directive']
+
+The IPython directive outputs code-blocks with the language 'ipython'. So
+if you do not have the syntax highlighting extension enabled as well, then
+all rendered code-blocks will be uncolored. By default this directive assumes
+that your prompts are unchanged IPython ones, but this can be customized.
+The configurable options that can be placed in conf.py are:
+
+ipython_savefig_dir:
+ The directory in which to save the figures. This is relative to the
+ Sphinx source directory. The default is `html_static_path`.
+ipython_rgxin:
+ The compiled regular expression to denote the start of IPython input
+ lines. The default is re.compile('In \[(\d+)\]:\s?(.*)\s*'). You
+ shouldn't need to change this.
+ipython_rgxout:
+ The compiled regular expression to denote the start of IPython output
+ lines. The default is re.compile('Out\[(\d+)\]:\s?(.*)\s*'). You
+ shouldn't need to change this.
+ipython_promptin:
+ The string to represent the IPython input prompt in the generated ReST.
+ The default is 'In [%d]:'. This expects that the line numbers are used
+ in the prompt.
+ipython_promptout:
+ The string to represent the IPython prompt in the generated ReST. The
+ default is 'Out [%d]:'. This expects that the line numbers are used
+ in the prompt.
+ipython_mplbackend:
+ The string which specifies if the embedded Sphinx shell should import
+ Matplotlib and set the backend. The value specifies a backend that is
+ passed to `matplotlib.use()` before any lines in `ipython_execlines` are
+ executed. If not specified in conf.py, then the default value of 'agg' is
+ used. To use the IPython directive without matplotlib as a dependency, set
+ the value to `None`. It may end up that matplotlib is still imported
+ if the user specifies so in `ipython_execlines` or makes use of the
+ @savefig pseudo decorator.
+ipython_execlines:
+ A list of strings to be exec'd in the embedded Sphinx shell. Typical
+ usage is to make certain packages always available. Set this to an empty
+ list if you wish to have no imports always available. If specified in
+ conf.py as `None`, then it has the effect of making no imports available.
+ If omitted from conf.py altogether, then the default value of
+ ['import numpy as np', 'import matplotlib.pyplot as plt'] is used.
+ipython_holdcount
+ When the @suppress pseudo-decorator is used, the execution count can be
+ incremented or not. The default behavior is to hold the execution count,
+ corresponding to a value of `True`. Set this to `False` to increment
+ the execution count after each suppressed command.
+
+As an example, to use the IPython directive when `matplotlib` is not available,
+one sets the backend to `None`::
+
+ ipython_mplbackend = None
+
+An example usage of the directive is:
+
+.. code-block:: rst
+
+ .. ipython::
+
+ In [1]: x = 1
+
+ In [2]: y = x**2
+
+ In [3]: print(y)
+
+See http://matplotlib.org/sampledoc/ipython_directive.html for additional
+documentation.
+
+Pseudo-Decorators
+=================
+
+Note: Only one decorator is supported per input. If more than one decorator
+is specified, then only the last one is used.
+
+In addition to the Pseudo-Decorators/options described at the above link,
+several enhancements have been made. The directive will emit a message to the
+console at build-time if code-execution resulted in an exception or warning.
+You can suppress these on a per-block basis by specifying the :okexcept:
+or :okwarning: options:
+
+.. code-block:: rst
+
+ .. ipython::
+ :okexcept:
+ :okwarning:
+
+ In [1]: 1/0
+ In [2]: # raise warning.
+
+ToDo
+----
+
+- Turn the ad-hoc test() function into a real test suite.
+- Break up ipython-specific functionality from matplotlib stuff into better
+ separated code.
+
+Authors
+-------
+
+- John D Hunter: orignal author.
+- Fernando Perez: refactoring, documentation, cleanups, port to 0.11.
+- VáclavŠmilauer <eudoxos-AT-arcig.cz>: Prompt generalizations.
+- Skipper Seabold, refactoring, cleanups, pure python addition
+"""
+from __future__ import print_function
+
+#-----------------------------------------------------------------------------
+# Imports
+#-----------------------------------------------------------------------------
+
+# Stdlib
+import os
+import re
+import sys
+import tempfile
+import ast
+import warnings
+
+# To keep compatibility with various python versions
+try:
+ from hashlib import md5
+except ImportError:
+ from md5 import md5
+
+# Third-party
+import sphinx
+from docutils.parsers.rst import directives
+from docutils import nodes
+from sphinx.util.compat import Directive
+
+# Our own
+from IPython import Config, InteractiveShell
+from IPython.core.profiledir import ProfileDir
+from IPython.utils import io
+from IPython.utils.py3compat import PY3
+
+if PY3:
+ from io import StringIO
+else:
+ from StringIO import StringIO
+
+#-----------------------------------------------------------------------------
+# Globals
+#-----------------------------------------------------------------------------
+# for tokenizing blocks
+COMMENT, INPUT, OUTPUT = range(3)
+
+#-----------------------------------------------------------------------------
+# Functions and class declarations
+#-----------------------------------------------------------------------------
+
+def block_parser(part, rgxin, rgxout, fmtin, fmtout):
+ """
+ part is a string of ipython text, comprised of at most one
+ input, one ouput, comments, and blank lines. The block parser
+ parses the text into a list of::
+
+ blocks = [ (TOKEN0, data0), (TOKEN1, data1), ...]
+
+ where TOKEN is one of [COMMENT | INPUT | OUTPUT ] and
+ data is, depending on the type of token::
+
+ COMMENT : the comment string
+
+ INPUT: the (DECORATOR, INPUT_LINE, REST) where
+ DECORATOR: the input decorator (or None)
+ INPUT_LINE: the input as string (possibly multi-line)
+ REST : any stdout generated by the input line (not OUTPUT)
+
+ OUTPUT: the output string, possibly multi-line
+
+ """
+ block = []
+ lines = part.split('\n')
+ N = len(lines)
+ i = 0
+ decorator = None
+ while 1:
+
+ if i==N:
+ # nothing left to parse -- the last line
+ break
+
+ line = lines[i]
+ i += 1
+ line_stripped = line.strip()
+ if line_stripped.startswith('#'):
+ block.append((COMMENT, line))
+ continue
+
+ if line_stripped.startswith('@'):
+ # Here is where we assume there is, at most, one decorator.
+ # Might need to rethink this.
+ decorator = line_stripped
+ continue
+
+ # does this look like an input line?
+ matchin = rgxin.match(line)
+ if matchin:
+ lineno, inputline = int(matchin.group(1)), matchin.group(2)
+
+ # the ....: continuation string
+ continuation = ' %s:'%''.join(['.']*(len(str(lineno))+2))
+ Nc = len(continuation)
+ # input lines can continue on for more than one line, if
+ # we have a '\' line continuation char or a function call
+ # echo line 'print'. The input line can only be
+ # terminated by the end of the block or an output line, so
+ # we parse out the rest of the input line if it is
+ # multiline as well as any echo text
+
+ rest = []
+ while i<N:
+
+ # look ahead; if the next line is blank, or a comment, or
+ # an output line, we're done
+
+ nextline = lines[i]
+ matchout = rgxout.match(nextline)
+ #print "nextline=%s, continuation=%s, starts=%s"%(nextline, continuation, nextline.startswith(continuation))
+ if matchout or nextline.startswith('#'):
+ break
+ elif nextline.startswith(continuation):
+ # The default ipython_rgx* treat the space following the colon as optional.
+ # However, If the space is there we must consume it or code
+ # employing the cython_magic extension will fail to execute.
+ #
+ # This works with the default ipython_rgx* patterns,
+ # If you modify them, YMMV.
+ nextline = nextline[Nc:]
+ if nextline and nextline[0] == ' ':
+ nextline = nextline[1:]
+
+ inputline += '\n' + nextline
+ else:
+ rest.append(nextline)
+ i+= 1
+
+ block.append((INPUT, (decorator, inputline, '\n'.join(rest))))
+ continue
+
+ # if it looks like an output line grab all the text to the end
+ # of the block
+ matchout = rgxout.match(line)
+ if matchout:
+ lineno, output = int(matchout.group(1)), matchout.group(2)
+ if i<N-1:
+ output = '\n'.join([output] + lines[i:])
+
+ block.append((OUTPUT, output))
+ break
+
+ return block
+
+
+class EmbeddedSphinxShell(object):
+ """An embedded IPython instance to run inside Sphinx"""
+
+ def __init__(self, exec_lines=None):
+
+ self.cout = StringIO()
+
+ if exec_lines is None:
+ exec_lines = []
+
+ # Create config object for IPython
+ config = Config()
+ config.InteractiveShell.autocall = False
+ config.InteractiveShell.autoindent = False
+ config.InteractiveShell.colors = 'NoColor'
+
+ # create a profile so instance history isn't saved
+ tmp_profile_dir = tempfile.mkdtemp(prefix='profile_')
+ profname = 'auto_profile_sphinx_build'
+ pdir = os.path.join(tmp_profile_dir,profname)
+ profile = ProfileDir.create_profile_dir(pdir)
+
+ # Create and initialize global ipython, but don't start its mainloop.
+ # This will persist across different EmbededSphinxShell instances.
+ IP = InteractiveShell.instance(config=config, profile_dir=profile)
+
+ # io.stdout redirect must be done after instantiating InteractiveShell
+ io.stdout = self.cout
+ io.stderr = self.cout
+
+ # For debugging, so we can see normal output, use this:
+ #from IPython.utils.io import Tee
+ #io.stdout = Tee(self.cout, channel='stdout') # dbg
+ #io.stderr = Tee(self.cout, channel='stderr') # dbg
+
+ # Store a few parts of IPython we'll need.
+ self.IP = IP
+ self.user_ns = self.IP.user_ns
+ self.user_global_ns = self.IP.user_global_ns
+
+ self.input = ''
+ self.output = ''
+
+ self.is_verbatim = False
+ self.is_doctest = False
+ self.is_suppress = False
+
+ # Optionally, provide more detailed information to shell.
+ # this is assigned by the SetUp method of IPythonDirective
+ # to point at itself.
+ #
+ # So, you can access handy things at self.directive.state
+ self.directive = None
+
+ # on the first call to the savefig decorator, we'll import
+ # pyplot as plt so we can make a call to the plt.gcf().savefig
+ self._pyplot_imported = False
+
+ # Prepopulate the namespace.
+ for line in exec_lines:
+ self.process_input_line(line, store_history=False)
+
+ def clear_cout(self):
+ self.cout.seek(0)
+ self.cout.truncate(0)
+
+ def process_input_line(self, line, store_history=True):
+ """process the input, capturing stdout"""
+
+ stdout = sys.stdout
+ splitter = self.IP.input_splitter
+ try:
+ sys.stdout = self.cout
+ splitter.push(line)
+ more = splitter.push_accepts_more()
+ if not more:
+ source_raw = splitter.raw_reset()
+ self.IP.run_cell(source_raw, store_history=store_history)
+ finally:
+ sys.stdout = stdout
+
+ def process_image(self, decorator):
+ """
+ # build out an image directive like
+ # .. image:: somefile.png
+ # :width 4in
+ #
+ # from an input like
+ # savefig somefile.png width=4in
+ """
+ savefig_dir = self.savefig_dir
+ source_dir = self.source_dir
+ saveargs = decorator.split(' ')
+ filename = saveargs[1]
+ # insert relative path to image file in source
+ outfile = os.path.relpath(os.path.join(savefig_dir,filename),
+ source_dir)
+
+ imagerows = ['.. image:: %s'%outfile]
+
+ for kwarg in saveargs[2:]:
+ arg, val = kwarg.split('=')
+ arg = arg.strip()
+ val = val.strip()
+ imagerows.append(' :%s: %s'%(arg, val))
+
+ image_file = os.path.basename(outfile) # only return file name
+ image_directive = '\n'.join(imagerows)
+ return image_file, image_directive
+
+ # Callbacks for each type of token
+ def process_input(self, data, input_prompt, lineno):
+ """
+ Process data block for INPUT token.
+
+ """
+ decorator, input, rest = data
+ image_file = None
+ image_directive = None
+
+ is_verbatim = decorator=='@verbatim' or self.is_verbatim
+ is_doctest = (decorator is not None and \
+ decorator.startswith('@doctest')) or self.is_doctest
+ is_suppress = decorator=='@suppress' or self.is_suppress
+ is_okexcept = decorator=='@okexcept' or self.is_okexcept
+ is_okwarning = decorator=='@okwarning' or self.is_okwarning
+ is_savefig = decorator is not None and \
+ decorator.startswith('@savefig')
+
+ input_lines = input.split('\n')
+ if len(input_lines) > 1:
+ if input_lines[-1] != "":
+ input_lines.append('') # make sure there's a blank line
+ # so splitter buffer gets reset
+
+ continuation = ' %s:'%''.join(['.']*(len(str(lineno))+2))
+
+ if is_savefig:
+ image_file, image_directive = self.process_image(decorator)
+
+ ret = []
+ is_semicolon = False
+
+ # Hold the execution count, if requested to do so.
+ if is_suppress and self.hold_count:
+ store_history = False
+ else:
+ store_history = True
+
+ # Note: catch_warnings is not thread safe
+ with warnings.catch_warnings(record=True) as ws:
+ for i, line in enumerate(input_lines):
+ if line.endswith(';'):
+ is_semicolon = True
+
+ if i == 0:
+ # process the first input line
+ if is_verbatim:
+ self.process_input_line('')
+ self.IP.execution_count += 1 # increment it anyway
+ else:
+ # only submit the line in non-verbatim mode
+ self.process_input_line(line, store_history=store_history)
+ formatted_line = '%s %s'%(input_prompt, line)
+ else:
+ # process a continuation line
+ if not is_verbatim:
+ self.process_input_line(line, store_history=store_history)
+
+ formatted_line = '%s %s'%(continuation, line)
+
+ if not is_suppress:
+ ret.append(formatted_line)
+
+ if not is_suppress and len(rest.strip()) and is_verbatim:
+ # The "rest" is the standard output of the input. This needs to be
+ # added when in verbatim mode. If there is no "rest", then we don't
+ # add it, as the new line will be added by the processed output.
+ ret.append(rest)
+
+ # Fetch the processed output. (This is not the submitted output.)
+ self.cout.seek(0)
+ processed_output = self.cout.read()
+ if not is_suppress and not is_semicolon:
+ #
+ # In IPythonDirective.run, the elements of `ret` are eventually
+ # combined such that '' entries correspond to newlines. So if
+ # `processed_output` is equal to '', then the adding it to `ret`
+ # ensures that there is a blank line between consecutive inputs
+ # that have no outputs, as in:
+ #
+ # In [1]: x = 4
+ #
+ # In [2]: x = 5
+ #
+ # When there is processed output, it has a '\n' at the tail end. So
+ # adding the output to `ret` will provide the necessary spacing
+ # between consecutive input/output blocks, as in:
+ #
+ # In [1]: x
+ # Out[1]: 5
+ #
+ # In [2]: x
+ # Out[2]: 5
+ #
+ # When there is stdout from the input, it also has a '\n' at the
+ # tail end, and so this ensures proper spacing as well. E.g.:
+ #
+ # In [1]: print x
+ # 5
+ #
+ # In [2]: x = 5
+ #
+ # When in verbatim mode, `processed_output` is empty (because
+ # nothing was passed to IP. Sometimes the submitted code block has
+ # an Out[] portion and sometimes it does not. When it does not, we
+ # need to ensure proper spacing, so we have to add '' to `ret`.
+ # However, if there is an Out[] in the submitted code, then we do
+ # not want to add a newline as `process_output` has stuff to add.
+ # The difficulty is that `process_input` doesn't know if
+ # `process_output` will be called---so it doesn't know if there is
+ # Out[] in the code block. The requires that we include a hack in
+ # `process_block`. See the comments there.
+ #
+ ret.append(processed_output)
+ elif is_semicolon:
+ # Make sure there is a newline after the semicolon.
+ ret.append('')
+
+ # context information
+ filename = "Unknown"
+ lineno = 0
+ if self.directive.state:
+ filename = self.directive.state.document.current_source
+ lineno = self.directive.state.document.current_line
+
+ # output any exceptions raised during execution to stdout
+ # unless :okexcept: has been specified.
+ if not is_okexcept and "Traceback" in processed_output:
+ s = "\nException in %s at block ending on line %s\n" % (filename, lineno)
+ s += "Specify :okexcept: as an option in the ipython:: block to suppress this message\n"
+ sys.stdout.write('\n\n>>>' + ('-' * 73))
+ sys.stdout.write(s)
+ sys.stdout.write(processed_output)
+ sys.stdout.write('<<<' + ('-' * 73) + '\n\n')
+
+ # output any warning raised during execution to stdout
+ # unless :okwarning: has been specified.
+ if not is_okwarning:
+ for w in ws:
+ s = "\nWarning in %s at block ending on line %s\n" % (filename, lineno)
+ s += "Specify :okwarning: as an option in the ipython:: block to suppress this message\n"
+ sys.stdout.write('\n\n>>>' + ('-' * 73))
+ sys.stdout.write(s)
+ sys.stdout.write(('-' * 76) + '\n')
+ s=warnings.formatwarning(w.message, w.category,
+ w.filename, w.lineno, w.line)
+ sys.stdout.write(s)
+ sys.stdout.write('<<<' + ('-' * 73) + '\n')
+
+ self.cout.truncate(0)
+
+ return (ret, input_lines, processed_output,
+ is_doctest, decorator, image_file, image_directive)
+
+
+ def process_output(self, data, output_prompt, input_lines, output,
+ is_doctest, decorator, image_file):
+ """
+ Process data block for OUTPUT token.
+
+ """
+ # Recall: `data` is the submitted output, and `output` is the processed
+ # output from `input_lines`.
+
+ TAB = ' ' * 4
+
+ if is_doctest and output is not None:
+
+ found = output # This is the processed output
+ found = found.strip()
+ submitted = data.strip()
+
+ if self.directive is None:
+ source = 'Unavailable'
+ content = 'Unavailable'
+ else:
+ source = self.directive.state.document.current_source
+ content = self.directive.content
+ # Add tabs and join into a single string.
+ content = '\n'.join([TAB + line for line in content])
+
+ # Make sure the output contains the output prompt.
+ ind = found.find(output_prompt)
+ if ind < 0:
+ e = ('output does not contain output prompt\n\n'
+ 'Document source: {0}\n\n'
+ 'Raw content: \n{1}\n\n'
+ 'Input line(s):\n{TAB}{2}\n\n'
+ 'Output line(s):\n{TAB}{3}\n\n')
+ e = e.format(source, content, '\n'.join(input_lines),
+ repr(found), TAB=TAB)
+ raise RuntimeError(e)
+ found = found[len(output_prompt):].strip()
+
+ # Handle the actual doctest comparison.
+ if decorator.strip() == '@doctest':
+ # Standard doctest
+ if found != submitted:
+ e = ('doctest failure\n\n'
+ 'Document source: {0}\n\n'
+ 'Raw content: \n{1}\n\n'
+ 'On input line(s):\n{TAB}{2}\n\n'
+ 'we found output:\n{TAB}{3}\n\n'
+ 'instead of the expected:\n{TAB}{4}\n\n')
+ e = e.format(source, content, '\n'.join(input_lines),
+ repr(found), repr(submitted), TAB=TAB)
+ raise RuntimeError(e)
+ else:
+ self.custom_doctest(decorator, input_lines, found, submitted)
+
+ # When in verbatim mode, this holds additional submitted output
+ # to be written in the final Sphinx output.
+ # https://github.com/ipython/ipython/issues/5776
+ out_data = []
+
+ is_verbatim = decorator=='@verbatim' or self.is_verbatim
+ if is_verbatim and data.strip():
+ # Note that `ret` in `process_block` has '' as its last element if
+ # the code block was in verbatim mode. So if there is no submitted
+ # output, then we will have proper spacing only if we do not add
+ # an additional '' to `out_data`. This is why we condition on
+ # `and data.strip()`.
+
+ # The submitted output has no output prompt. If we want the
+ # prompt and the code to appear, we need to join them now
+ # instead of adding them separately---as this would create an
+ # undesired newline. How we do this ultimately depends on the
+ # format of the output regex. I'll do what works for the default
+ # prompt for now, and we might have to adjust if it doesn't work
+ # in other cases. Finally, the submitted output does not have
+ # a trailing newline, so we must add it manually.
+ out_data.append("{0} {1}\n".format(output_prompt, data))
+
+ return out_data
+
+ def process_comment(self, data):
+ """Process data fPblock for COMMENT token."""
+ if not self.is_suppress:
+ return [data]
+
+ def save_image(self, image_file):
+ """
+ Saves the image file to disk.
+ """
+ self.ensure_pyplot()
+ command = 'plt.gcf().savefig("%s")'%image_file
+ #print 'SAVEFIG', command # dbg
+ self.process_input_line('bookmark ipy_thisdir', store_history=False)
+ self.process_input_line('cd -b ipy_savedir', store_history=False)
+ self.process_input_line(command, store_history=False)
+ self.process_input_line('cd -b ipy_thisdir', store_history=False)
+ self.process_input_line('bookmark -d ipy_thisdir', store_history=False)
+ self.clear_cout()
+
+ def process_block(self, block):
+ """
+ process block from the block_parser and return a list of processed lines
+ """
+ ret = []
+ output = None
+ input_lines = None
+ lineno = self.IP.execution_count
+
+ input_prompt = self.promptin % lineno
+ output_prompt = self.promptout % lineno
+ image_file = None
+ image_directive = None
+
+ for token, data in block:
+ if token == COMMENT:
+ out_data = self.process_comment(data)
+ elif token == INPUT:
+ (out_data, input_lines, output, is_doctest,
+ decorator, image_file, image_directive) = \
+ self.process_input(data, input_prompt, lineno)
+ elif token == OUTPUT:
+ out_data = \
+ self.process_output(data, output_prompt, input_lines,
+ output, is_doctest, decorator,
+ image_file)
+ if out_data:
+ # Then there was user submitted output in verbatim mode.
+ # We need to remove the last element of `ret` that was
+ # added in `process_input`, as it is '' and would introduce
+ # an undesirable newline.
+ assert(ret[-1] == '')
+ del ret[-1]
+
+ if out_data:
+ ret.extend(out_data)
+
+ # save the image files
+ if image_file is not None:
+ self.save_image(image_file)
+
+ return ret, image_directive
+
+ def ensure_pyplot(self):
+ """
+ Ensures that pyplot has been imported into the embedded IPython shell.
+
+ Also, makes sure to set the backend appropriately if not set already.
+
+ """
+ # We are here if the @figure pseudo decorator was used. Thus, it's
+ # possible that we could be here even if python_mplbackend were set to
+ # `None`. That's also strange and perhaps worthy of raising an
+ # exception, but for now, we just set the backend to 'agg'.
+
+ if not self._pyplot_imported:
+ if 'matplotlib.backends' not in sys.modules:
+ # Then ipython_matplotlib was set to None but there was a
+ # call to the @figure decorator (and ipython_execlines did
+ # not set a backend).
+ #raise Exception("No backend was set, but @figure was used!")
+ import matplotlib
+ matplotlib.use('agg')
+
+ # Always import pyplot into embedded shell.
+ self.process_input_line('import matplotlib.pyplot as plt',
+ store_history=False)
+ self._pyplot_imported = True
+
+ def process_pure_python(self, content):
+ """
+ content is a list of strings. it is unedited directive content
+
+ This runs it line by line in the InteractiveShell, prepends
+ prompts as needed capturing stderr and stdout, then returns
+ the content as a list as if it were ipython code
+ """
+ output = []
+ savefig = False # keep up with this to clear figure
+ multiline = False # to handle line continuation
+ multiline_start = None
+ fmtin = self.promptin
+
+ ct = 0
+
+ for lineno, line in enumerate(content):
+
+ line_stripped = line.strip()
+ if not len(line):
+ output.append(line)
+ continue
+
+ # handle decorators
+ if line_stripped.startswith('@'):
+ output.extend([line])
+ if 'savefig' in line:
+ savefig = True # and need to clear figure
+ continue
+
+ # handle comments
+ if line_stripped.startswith('#'):
+ output.extend([line])
+ continue
+
+ # deal with lines checking for multiline
+ continuation = u' %s:'% ''.join(['.']*(len(str(ct))+2))
+ if not multiline:
+ modified = u"%s %s" % (fmtin % ct, line_stripped)
+ output.append(modified)
+ ct += 1
+ try:
+ ast.parse(line_stripped)
+ output.append(u'')
+ except Exception: # on a multiline
+ multiline = True
+ multiline_start = lineno
+ else: # still on a multiline
+ modified = u'%s %s' % (continuation, line)
+ output.append(modified)
+
+ # if the next line is indented, it should be part of multiline
+ if len(content) > lineno + 1:
+ nextline = content[lineno + 1]
+ if len(nextline) - len(nextline.lstrip()) > 3:
+ continue
+ try:
+ mod = ast.parse(
+ '\n'.join(content[multiline_start:lineno+1]))
+ if isinstance(mod.body[0], ast.FunctionDef):
+ # check to see if we have the whole function
+ for element in mod.body[0].body:
+ if isinstance(element, ast.Return):
+ multiline = False
+ else:
+ output.append(u'')
+ multiline = False
+ except Exception:
+ pass
+
+ if savefig: # clear figure if plotted
+ self.ensure_pyplot()
+ self.process_input_line('plt.clf()', store_history=False)
+ self.clear_cout()
+ savefig = False
+
+ return output
+
+ def custom_doctest(self, decorator, input_lines, found, submitted):
+ """
+ Perform a specialized doctest.
+
+ """
+ from .custom_doctests import doctests
+
+ args = decorator.split()
+ doctest_type = args[1]
+ if doctest_type in doctests:
+ doctests[doctest_type](self, args, input_lines, found, submitted)
+ else:
+ e = "Invalid option to @doctest: {0}".format(doctest_type)
+ raise Exception(e)
+
+
+class IPythonDirective(Directive):
+
+ has_content = True
+ required_arguments = 0
+ optional_arguments = 4 # python, suppress, verbatim, doctest
+ final_argumuent_whitespace = True
+ option_spec = { 'python': directives.unchanged,
+ 'suppress' : directives.flag,
+ 'verbatim' : directives.flag,
+ 'doctest' : directives.flag,
+ 'okexcept': directives.flag,
+ 'okwarning': directives.flag
+ }
+
+ shell = None
+
+ seen_docs = set()
+
+ def get_config_options(self):
+ # contains sphinx configuration variables
+ config = self.state.document.settings.env.config
+
+ # get config variables to set figure output directory
+ confdir = self.state.document.settings.env.app.confdir
+ savefig_dir = config.ipython_savefig_dir
+ source_dir = os.path.dirname(self.state.document.current_source)
+ if savefig_dir is None:
+ savefig_dir = config.html_static_path
+ if isinstance(savefig_dir, list):
+ savefig_dir = savefig_dir[0] # safe to assume only one path?
+ savefig_dir = os.path.join(confdir, savefig_dir)
+
+ # get regex and prompt stuff
+ rgxin = config.ipython_rgxin
+ rgxout = config.ipython_rgxout
+ promptin = config.ipython_promptin
+ promptout = config.ipython_promptout
+ mplbackend = config.ipython_mplbackend
+ exec_lines = config.ipython_execlines
+ hold_count = config.ipython_holdcount
+
+ return (savefig_dir, source_dir, rgxin, rgxout,
+ promptin, promptout, mplbackend, exec_lines, hold_count)
+
+ def setup(self):
+ # Get configuration values.
+ (savefig_dir, source_dir, rgxin, rgxout, promptin, promptout,
+ mplbackend, exec_lines, hold_count) = self.get_config_options()
+
+ if self.shell is None:
+ # We will be here many times. However, when the
+ # EmbeddedSphinxShell is created, its interactive shell member
+ # is the same for each instance.
+
+ if mplbackend:
+ import matplotlib
+ # Repeated calls to use() will not hurt us since `mplbackend`
+ # is the same each time.
+ matplotlib.use(mplbackend)
+
+ # Must be called after (potentially) importing matplotlib and
+ # setting its backend since exec_lines might import pylab.
+ self.shell = EmbeddedSphinxShell(exec_lines)
+
+ # Store IPython directive to enable better error messages
+ self.shell.directive = self
+
+ # reset the execution count if we haven't processed this doc
+ #NOTE: this may be borked if there are multiple seen_doc tmp files
+ #check time stamp?
+ if not self.state.document.current_source in self.seen_docs:
+ self.shell.IP.history_manager.reset()
+ self.shell.IP.execution_count = 1
+ self.shell.IP.prompt_manager.width = 0
+ self.seen_docs.add(self.state.document.current_source)
+
+ # and attach to shell so we don't have to pass them around
+ self.shell.rgxin = rgxin
+ self.shell.rgxout = rgxout
+ self.shell.promptin = promptin
+ self.shell.promptout = promptout
+ self.shell.savefig_dir = savefig_dir
+ self.shell.source_dir = source_dir
+ self.shell.hold_count = hold_count
+
+ # setup bookmark for saving figures directory
+ self.shell.process_input_line('bookmark ipy_savedir %s'%savefig_dir,
+ store_history=False)
+ self.shell.clear_cout()
+
+ return rgxin, rgxout, promptin, promptout
+
+ def teardown(self):
+ # delete last bookmark
+ self.shell.process_input_line('bookmark -d ipy_savedir',
+ store_history=False)
+ self.shell.clear_cout()
+
+ def run(self):
+ debug = False
+
+ #TODO, any reason block_parser can't be a method of embeddable shell
+ # then we wouldn't have to carry these around
+ rgxin, rgxout, promptin, promptout = self.setup()
+
+ options = self.options
+ self.shell.is_suppress = 'suppress' in options
+ self.shell.is_doctest = 'doctest' in options
+ self.shell.is_verbatim = 'verbatim' in options
+ self.shell.is_okexcept = 'okexcept' in options
+ self.shell.is_okwarning = 'okwarning' in options
+
+ # handle pure python code
+ if 'python' in self.arguments:
+ content = self.content
+ self.content = self.shell.process_pure_python(content)
+
+ parts = '\n'.join(self.content).split('\n\n')
+
+ lines = ['.. code-block:: ipython', '']
+ figures = []
+
+ for part in parts:
+ block = block_parser(part, rgxin, rgxout, promptin, promptout)
+ if len(block):
+ rows, figure = self.shell.process_block(block)
+ for row in rows:
+ lines.extend([' {0}'.format(line)
+ for line in row.split('\n')])
+
+ if figure is not None:
+ figures.append(figure)
+
+ for figure in figures:
+ lines.append('')
+ lines.extend(figure.split('\n'))
+ lines.append('')
+
+ if len(lines) > 2:
+ if debug:
+ print('\n'.join(lines))
+ else:
+ # This has to do with input, not output. But if we comment
+ # these lines out, then no IPython code will appear in the
+ # final output.
+ self.state_machine.insert_input(
+ lines, self.state_machine.input_lines.source(0))
+
+ # cleanup
+ self.teardown()
+
+ return []
+
+# Enable as a proper Sphinx directive
+def setup(app):
+ setup.app = app
+
+ app.add_directive('ipython', IPythonDirective)
+ app.add_config_value('ipython_savefig_dir', None, 'env')
+ app.add_config_value('ipython_rgxin',
+ re.compile('In \[(\d+)\]:\s?(.*)\s*'), 'env')
+ app.add_config_value('ipython_rgxout',
+ re.compile('Out\[(\d+)\]:\s?(.*)\s*'), 'env')
+ app.add_config_value('ipython_promptin', 'In [%d]:', 'env')
+ app.add_config_value('ipython_promptout', 'Out[%d]:', 'env')
+
+ # We could just let matplotlib pick whatever is specified as the default
+ # backend in the matplotlibrc file, but this would cause issues if the
+ # backend didn't work in headless environments. For this reason, 'agg'
+ # is a good default backend choice.
+ app.add_config_value('ipython_mplbackend', 'agg', 'env')
+
+ # If the user sets this config value to `None`, then EmbeddedSphinxShell's
+ # __init__ method will treat it as [].
+ execlines = ['import numpy as np', 'import matplotlib.pyplot as plt']
+ app.add_config_value('ipython_execlines', execlines, 'env')
+
+ app.add_config_value('ipython_holdcount', True, 'env')
+
+# Simple smoke test, needs to be converted to a proper automatic test.
+def test():
+
+ examples = [
+ r"""
+In [9]: pwd
+Out[9]: '/home/jdhunter/py4science/book'
+
+In [10]: cd bookdata/
+/home/jdhunter/py4science/book/bookdata
+
+In [2]: from pylab import *
+
+In [2]: ion()
+
+In [3]: im = imread('stinkbug.png')
+
+@savefig mystinkbug.png width=4in
+In [4]: imshow(im)
+Out[4]: <matplotlib.image.AxesImage object at 0x39ea850>
+
+""",
+ r"""
+
+In [1]: x = 'hello world'
+
+# string methods can be
+# used to alter the string
+@doctest
+In [2]: x.upper()
+Out[2]: 'HELLO WORLD'
+
+@verbatim
+In [3]: x.st<TAB>
+x.startswith x.strip
+""",
+ r"""
+
+In [130]: url = 'http://ichart.finance.yahoo.com/table.csv?s=CROX\
+ .....: &d=9&e=22&f=2009&g=d&a=1&br=8&c=2006&ignore=.csv'
+
+In [131]: print url.split('&')
+['http://ichart.finance.yahoo.com/table.csv?s=CROX', 'd=9', 'e=22', 'f=2009', 'g=d', 'a=1', 'b=8', 'c=2006', 'ignore=.csv']
+
+In [60]: import urllib
+
+""",
+ r"""\
+
+In [133]: import numpy.random
+
+@suppress
+In [134]: numpy.random.seed(2358)
+
+@doctest
+In [135]: numpy.random.rand(10,2)
+Out[135]:
+array([[ 0.64524308, 0.59943846],
+ [ 0.47102322, 0.8715456 ],
+ [ 0.29370834, 0.74776844],
+ [ 0.99539577, 0.1313423 ],
+ [ 0.16250302, 0.21103583],
+ [ 0.81626524, 0.1312433 ],
+ [ 0.67338089, 0.72302393],
+ [ 0.7566368 , 0.07033696],
+ [ 0.22591016, 0.77731835],
+ [ 0.0072729 , 0.34273127]])
+
+""",
+
+ r"""
+In [106]: print x
+jdh
+
+In [109]: for i in range(10):
+ .....: print i
+ .....:
+ .....:
+0
+1
+2
+3
+4
+5
+6
+7
+8
+9
+""",
+
+ r"""
+
+In [144]: from pylab import *
+
+In [145]: ion()
+
+# use a semicolon to suppress the output
+@savefig test_hist.png width=4in
+In [151]: hist(np.random.randn(10000), 100);
+
+
+@savefig test_plot.png width=4in
+In [151]: plot(np.random.randn(10000), 'o');
+ """,
+
+ r"""
+# use a semicolon to suppress the output
+In [151]: plt.clf()
+
+@savefig plot_simple.png width=4in
+In [151]: plot([1,2,3])
+
+@savefig hist_simple.png width=4in
+In [151]: hist(np.random.randn(10000), 100);
+
+""",
+ r"""
+# update the current fig
+In [151]: ylabel('number')
+
+In [152]: title('normal distribution')
+
+
+@savefig hist_with_text.png
+In [153]: grid(True)
+
+@doctest float
+In [154]: 0.1 + 0.2
+Out[154]: 0.3
+
+@doctest float
+In [155]: np.arange(16).reshape(4,4)
+Out[155]:
+array([[ 0, 1, 2, 3],
+ [ 4, 5, 6, 7],
+ [ 8, 9, 10, 11],
+ [12, 13, 14, 15]])
+
+In [1]: x = np.arange(16, dtype=float).reshape(4,4)
+
+In [2]: x[0,0] = np.inf
+
+In [3]: x[0,1] = np.nan
+
+@doctest float
+In [4]: x
+Out[4]:
+array([[ inf, nan, 2., 3.],
+ [ 4., 5., 6., 7.],
+ [ 8., 9., 10., 11.],
+ [ 12., 13., 14., 15.]])
+
+
+ """,
+ ]
+ # skip local-file depending first example:
+ examples = examples[1:]
+
+ #ipython_directive.DEBUG = True # dbg
+ #options = dict(suppress=True) # dbg
+ options = dict()
+ for example in examples:
+ content = example.split('\n')
+ IPythonDirective('debug', arguments=None, options=options,
+ content=content, lineno=0,
+ content_offset=None, block_text=None,
+ state=None, state_machine=None,
+ )
+
+# Run test suite as a script
+if __name__=='__main__':
+ if not os.path.isdir('_static'):
+ os.mkdir('_static')
+ test()
+ print('All OK? Check figures in _static/')
=== modified file 'examples/capillary/liquidmigration/showcase.py'
--- examples/capillary/liquidmigration/showcase.py 2014-05-08 05:57:04 +0000
+++ examples/capillary/liquidmigration/showcase.py 2014-06-18 07:35:57 +0000
@@ -29,8 +29,8 @@
id5 = O.bodies.append(sphere(center=[(r1+r2)*d*2,(r1+r2)*d,0], radius=r2,material=mat1,fixed=True, color=[0,1,0]))
-Vf = 0.5e-1
-Vfmin = 0.1e-1
+Vf = 0.0e-1
+Vfmin = 0.0e-1
O.bodies[id1].Vf = Vf
O.bodies[id1].Vmin = Vfmin
@@ -61,7 +61,7 @@
[Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],
[Law2_ScGeom_ViscElCapPhys_Basic()],
),
- LiqControl(),
+ LiqControl(label='lqc'),
NewtonIntegrator(damping=0,gravity=[0,0,0]),
PyRunner(command='showData()',iterPeriod=1),
]
@@ -105,14 +105,18 @@
O.bodies[i].Vf = 0
O.bodies[i].Vmin = 0
-O.interactions[id1,id2].phys.Vb = 1.0
O.interactions[id1,id2].phys.Vmax = 5.0
-O.interactions[id2,id3].phys.Vb = 1.0
+lqc.addLiqInter(id1, id2, 1.0)
+
O.interactions[id2,id3].phys.Vmax = 5.0
-O.interactions[id3,id4].phys.Vb = 1.0
+lqc.addLiqInter(id2, id3, 1.0)
+
O.interactions[id3,id4].phys.Vmax = 5.0
-O.interactions[id3,id5].phys.Vb = 1.0
+lqc.addLiqInter(id3, id4, 1.0)
+
O.interactions[id3,id5].phys.Vmax = 5.0
+lqc.addLiqInter(id3, id5, 1.0)
+
O.run(1, True)
=== modified file 'examples/concrete/uniax.py'
--- examples/concrete/uniax.py 2013-08-09 02:27:56 +0000
+++ examples/concrete/uniax.py 2014-06-10 03:36:46 +0000
@@ -2,7 +2,7 @@
# -*- coding: utf-8 -*-
from __future__ import division
-from yade import plot,pack,timing,eudoxos
+from yade import plot,pack,timing
import time, sys, os, copy
#import matplotlib
@@ -142,8 +142,6 @@
extremum=max(sigma) if (strainer.strainRate>0) else min(sigma)
minMaxRatio=0.5 if mode=='tension' else 0.5
if extremum==0: return
- # uncomment to get graph for the very first time stopIfDamaged() is called
- #eudoxos.estimatePoissonYoung(principalAxis=axis,stress=strainer.avgStress,plot=True,cutoff=0.3)
import sys; sys.stdout.flush()
if abs(sigma[-1]/extremum)<minMaxRatio or abs(strainer.strain)>(5e-3 if isoPrestress==0 else 5e-2):
if mode=='tension' and doModes & 2: # only if compression is enabled
=== modified file 'examples/gts-horse/gts-random-pack-obb.py'
--- examples/gts-horse/gts-random-pack-obb.py 2013-08-09 02:27:56 +0000
+++ examples/gts-horse/gts-random-pack-obb.py 2014-06-12 16:54:31 +0000
@@ -13,7 +13,7 @@
poly=((3+.1,0),(3+0,.1),(3+sq2,.1+sq2),(3+.1+sq2,sq2),(3+.1,0))
#pylab.plot(*zip(*poly)); pylab.xlim(xmin=0); pylab.grid(); pylab.title('Meridian of the revolution surface\n(close to continue)'); pylab.gca().set_aspect(aspect='equal',adjustable='box'); pylab.show()
thetas=arange(0,pi/8,pi/24)
-pts=pack.revolutionSurfaceMeridians([poly for theta in thetas],thetas,origin=Vector3(-4,0,-1),orientation=Quaternion.Identity)
+pts=pack.revolutionSurfaceMeridians([poly for theta in thetas],thetas,origin=Vector3(-4,0,-1),orientation=Quaternion((0,1,0),0.0))
surf=pack.sweptPolylines2gtsSurface(pts,capStart=True,capEnd=True,threshold=1e-4)
O.bodies.append(pack.gtsSurface2Facets(surf,color=(1,0,1)))
# fill this solid with triaxial packing; it will compute minimum-volume oriented bounding box
=== modified file 'examples/test/performance/checkPerf.py'
--- examples/test/performance/checkPerf.py 2014-04-02 15:33:41 +0000
+++ examples/test/performance/checkPerf.py 2014-06-12 05:01:21 +0000
@@ -3,32 +3,40 @@
from yade import pack,export,geom,timing,bodiesHandling
import time,numpy
-radRAD=[23.658, #5000 elements
+radRAD=[
+ 23.658, #5000 elements
40.455, #25000 elements
50.97, #50000 elements
64.218, #100000 elements
- 80.91] #200000 elements
- #109.811] #500000 elements
+ 80.91, #200000 elements
+ #109.811, #500000 elements
+]
-iterN=[12000, #5000 elements
+iterN=[
+ 12000, #5000 elements
2500, #25000 elements
1400, #50000 elements
800, #100000 elements
- 200] #200000 elements
- #10] #500000 elements
+ 200, #200000 elements
+ #10, #500000 elements
+]
-coefCor=[110,
+coefCor=[
+ 110,
28,
18,
9,
- 2]
- #0.1]
+ 2,
+ #0.1,
+]
iterVel=[]
testTime=[]
particlesNumber=[]
+data=[]
numberTests = 3
+initIter = 1 # number of initial iterations excluded from performance check (e.g. Collider is initialised in first step therefore exclude first step from performance check)
tStartAll=time.time()
@@ -43,7 +51,6 @@
es=.003
frictionAngle=radians(35)
density=2300
-
defMat=O.materials.append(ViscElMat(density=density,frictionAngle=frictionAngle,tc=tc,en=en,et=es))
O.dt=.1*tc # time step
@@ -76,11 +83,11 @@
print "number of bodies %d"%len(O.bodies)
# Initialize the collider else it is not possible to compare the results with different nbIter
- O.run(1,1)
+ O.run(initIter,1)
O.timingEnabled=True
timing.reset()
tStart=time.time()
-
+ # run nbIter iterations
O.run(nbIter)
O.wait()
@@ -102,24 +109,56 @@
print
print
+print "___________________________________________________"
+print
+print "SUMMARY"
+print
scoreIterVel=0.0
for i in range(len(radRAD)):
iterAv=0.0
- iterVelNumpy=numpy.empty(3)
+ iterVelNumpy=numpy.empty(numberTests)
for z in range(numberTests):
iterVelNumpy[z]=iterVel[z*len(radRAD)+i]
avgVel = numpy.average(iterVelNumpy)
dispVel = numpy.std(iterVelNumpy)/numpy.average(iterVelNumpy)*100.0
- if (dispVel>10):
+ if (dispVel>10.):
print "Calculation velocity is unstable, try to close all programs and start performance tests again"
- print particlesNumber[i]," spheres, velocity=",avgVel, "+-",dispVel,"%"
+ print particlesNumber[i]," spheres, calculation velocity=",avgVel, "iter/sec +/-",dispVel,"%"
+ data+=[[particlesNumber[i],avgVel,avgVel*dispVel/100.]]
scoreIterVel+=avgVel/coefCor[i]*1000.0
print
-print
scoreIterVel = int(scoreIterVel)
-print "SCORE: " + str(scoreIterVel)
-print "Number of threads ", os.environ['OMP_NUM_THREADS']
-print"___________________________________________________"
-print "CPU info", os.system('cat /proc/cpuinfo')
+## write output file for graph
+import subprocess
+cmd = "cat /proc/cpuinfo | grep \'model name\' | uniq"
+#processor = subprocess.check_output(cmd, shell=True).lstrip('model name\t:').strip() # needs python >=2.7.0
+process = subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True)
+processor = process.communicate()[0].lstrip('model name\t:').strip()
+cmd = "cat /proc/cpuinfo | grep processor | wc -l"
+#cores = subprocess.check_output("cat /proc/cpuinfo | grep processor | wc -l", shell=True).strip() # needs python >=2.7.0
+process = subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True)
+cores = process.communicate()[0].strip()
+header='# '+ processor + ' ('+cores+' cores)'
+numThreads=os.environ['OMP_NUM_THREADS'] if (len(os.environ['OMP_NUM_THREADS'])==2) else ('0'+os.environ['OMP_NUM_THREADS'])
+filename=version+"_j"+numThreads+".dat"
+#numpy.savetxt(filename,numpy.array(data),fmt=['%i','%g','%g'], header=header) # needs numpy >=1.7.0
+fid = open( filename, 'w' )
+fid.write(header+"\n")
+for l in data:
+ fid.write("%d %f %f\n"%(l[0],l[1],l[2]))
+fid.close()
+print "Summary saved to "+filename
+print
+print "SCORE: " + str(scoreIterVel)
+print "Number of threads: ", os.environ['OMP_NUM_THREADS']
+print "___________________________________________________"
+print
+
+print "CPU info:"
+cmd = "lscpu"
+#cpuinfo=subprocess.check_output(cmd, shell=True) # needs python >=2.7.0
+process = subprocess.Popen(cmd, stdout=subprocess.PIPE, shell=True)
+cpuinfo = process.communicate()[0].lstrip('model name\t:').strip()
+print cpuinfo
sys.exit(0)
=== modified file 'gui/qt4/SerializableEditor.py'
--- gui/qt4/SerializableEditor.py 2013-03-18 20:03:32 +0000
+++ gui/qt4/SerializableEditor.py 2014-06-12 16:54:31 +0000
@@ -277,7 +277,7 @@
class Se3FakeType: pass
_fundamentalEditorMap={bool:AttrEditor_Bool,str:AttrEditor_Str,int:AttrEditor_Int,float:AttrEditor_Float,Quaternion:AttrEditor_Quaternion,Vector2:AttrEditor_Vector2,Vector3:AttrEditor_Vector3,Vector6:AttrEditor_Vector6,Matrix3:AttrEditor_Matrix3,Vector6i:AttrEditor_Vector6i,Vector3i:AttrEditor_Vector3i,Vector2i:AttrEditor_Vector2i,Se3FakeType:AttrEditor_Se3}
-_fundamentalInitValues={bool:True,str:'',int:0,float:0.0,Quaternion:Quaternion.Identity,Vector3:Vector3.Zero,Matrix3:Matrix3.Zero,Vector6:Vector6.Zero,Vector6i:Vector6i.Zero,Vector3i:Vector3i.Zero,Vector2i:Vector2i.Zero,Vector2:Vector2.Zero,Se3FakeType:(Vector3.Zero,Quaternion.Identity)}
+_fundamentalInitValues={bool:True,str:'',int:0,float:0.0,Quaternion:Quaternion((0,1,0),0.0),Vector3:Vector3.Zero,Matrix3:Matrix3.Zero,Vector6:Vector6.Zero,Vector6i:Vector6i.Zero,Vector3i:Vector3i.Zero,Vector2i:Vector2i.Zero,Vector2:Vector2.Zero,Se3FakeType:(Vector3.Zero,Quaternion((0,1,0),0.0))}
class SerQLabel(QLabel):
def __init__(self,parent,label,tooltip,path):
=== modified file 'lib/base/Math.cpp'
--- lib/base/Math.cpp 2010-11-07 11:46:20 +0000
+++ lib/base/Math.cpp 2014-06-17 10:18:00 +0000
@@ -10,3 +10,18 @@
template<> int ZeroInitializer<int>(){ return (int)0; }
template<> Real ZeroInitializer<Real>(){ return (Real)0; }
+
+#ifdef YADE_MASK_ARBITRARY
+bool operator==(const mask_t& g, int i) { return g == mask_t(i); }
+bool operator==(int i, const mask_t& g) { return g == i; }
+bool operator!=(const mask_t& g, int i) { return !(g == i); }
+bool operator!=(int i, const mask_t& g) { return g != i; }
+mask_t operator&(const mask_t& g, int i) { return g & mask_t(i); }
+mask_t operator&(int i, const mask_t& g) { return g & i; }
+mask_t operator|(const mask_t& g, int i) { return g | mask_t(i); }
+mask_t operator|(int i, const mask_t& g) { return g | i; }
+bool operator||(const mask_t& g, bool b) { return (g == 0) || b; }
+bool operator||(bool b, const mask_t& g) { return g || b; }
+bool operator&&(const mask_t& g, bool b) { return (g == 0) && b; }
+bool operator&&(bool b, const mask_t& g) { return g && b; }
+#endif
=== modified file 'lib/base/Math.hpp'
--- lib/base/Math.hpp 2013-08-23 15:21:20 +0000
+++ lib/base/Math.hpp 2014-06-17 10:18:00 +0000
@@ -8,16 +8,13 @@
typedef double Real;
#endif
+#define EIGEN_DONT_PARALLELIZE
+
#include<limits>
#include<cstdlib>
-
-// disable optimization which are "unsafe":
-// eigen objects cannot be passed by-value, otherwise they will no be aligned
-
-#define EIGEN_DONT_VECTORIZE
-#define EIGEN_DONT_ALIGN
-#define EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
-#define EIGEN_NO_DEBUG
+#ifdef YADE_MASK_ARBITRARY
+#include<bitset>
+#endif
#include<Eigen/Core>
#include<Eigen/Geometry>
@@ -209,12 +206,36 @@
template<typename Scalar> Scalar unitVectorsAngle(const VECTOR3_TEMPLATE(Scalar)& a, const VECTOR3_TEMPLATE(Scalar)& b){ return acos(a.dot(b)); }
// operators
+
+/*
+ * Mask
+ */
+#ifdef YADE_MASK_ARBITRARY
+typedef std::bitset<YADE_MASK_ARBITRARY_SIZE> mask_t;
+bool operator==(const mask_t& g, int i);
+bool operator==(int i, const mask_t& g);
+bool operator!=(const mask_t& g, int i);
+bool operator!=(int i, const mask_t& g);
+mask_t operator&(const mask_t& g, int i);
+mask_t operator&(int i, const mask_t& g);
+mask_t operator|(const mask_t& g, int i);
+mask_t operator|(int i, const mask_t& g);
+bool operator||(const mask_t& g, bool b);
+bool operator||(bool b, const mask_t& g);
+bool operator&&(const mask_t& g, bool b);
+bool operator&&(bool b, const mask_t& g);
+#else
+typedef int mask_t;
+#endif
+
+
/*
* typedefs
*/
typedef Se3<Real> Se3r;
+
/*
* Serialization of math classes
*/
@@ -313,6 +334,15 @@
BOOST_SERIALIZATION_NVP(m50) & BOOST_SERIALIZATION_NVP(m51) & BOOST_SERIALIZATION_NVP(m52) & BOOST_SERIALIZATION_NVP(m53) & BOOST_SERIALIZATION_NVP(m54) & BOOST_SERIALIZATION_NVP(m55);
}
+#ifdef YADE_MASK_ARBITRARY
+template<class Archive>
+void serialize(Archive & ar, mask_t& v, const unsigned int version){
+ std::string str = v.to_string();
+ ar & BOOST_SERIALIZATION_NVP(str);
+ v = mask_t(str);
+}
+#endif
+
} // namespace serialization
} // namespace boost
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2014-05-26 23:23:29 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2014-06-16 14:02:20 +0000
@@ -1068,7 +1068,7 @@
{
RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
static unsigned int number = 0;
- char filename[80];
+ char filename[250];
mkdir(folder, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
sprintf(filename,"%s/out_%d.vtk",folder,number++);
int firstReal=-1;
=== modified file 'pkg/common/InteractionLoop.cpp'
--- pkg/common/InteractionLoop.cpp 2014-05-23 13:05:19 +0000
+++ pkg/common/InteractionLoop.cpp 2014-06-06 19:26:47 +0000
@@ -70,7 +70,9 @@
const shared_ptr<Body>& b2_=Body::byId(I->getId2(),scene);
if(!b1_ || !b2_){ LOG_DEBUG("Body #"<<(b1_?I->getId2():I->getId1())<<" vanished, erasing intr #"<<I->getId1()<<"+#"<<I->getId2()<<"!"); scene->interactions->requestErase(I); continue; }
-
+
+ // Skip interaction with clumps
+ if (b1_->isClump() || b2_->isClump()) { continue; }
// we know there is no geometry functor already, take the short path
if(!I->functorCache.geomExists) { assert(!I->isReal()); continue; }
// no interaction geometry for either of bodies; no interaction possible
@@ -107,7 +109,7 @@
if(wasReal) scene->interactions->requestErase(I); // fully created interaction without geometry is reset and perhaps erased in the next step
continue; // in any case don't care about this one anymore
}
-
+
// IPhysDispatcher
if(!I->functorCache.phys){
I->functorCache.phys=physDispatcher->getFunctor2D(b1->material,b2->material,swap);
=== modified file 'pkg/dem/CapillaryTriaxialTest.cpp'
--- pkg/dem/CapillaryTriaxialTest.cpp 2014-06-05 13:19:44 +0000
+++ pkg/dem/CapillaryTriaxialTest.cpp 2014-06-17 10:18:00 +0000
@@ -209,7 +209,7 @@
void CapillaryTriaxialTest::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool big, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
+ body = shared_ptr<Body>(new Body); body->groupMask=2;
shared_ptr<FrictMat> physics(new FrictMat);
shared_ptr<Aabb> aabb(new Aabb);
@@ -254,7 +254,7 @@
void CapillaryTriaxialTest::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents, bool wire)
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
+ body = shared_ptr<Body>(new Body); body->groupMask=2;
shared_ptr<FrictMat> physics(new FrictMat);
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Box> iBox(new Box);
=== modified file 'pkg/dem/CohesiveTriaxialTest.cpp'
--- pkg/dem/CohesiveTriaxialTest.cpp 2014-06-05 13:19:44 +0000
+++ pkg/dem/CohesiveTriaxialTest.cpp 2014-06-17 10:18:00 +0000
@@ -200,7 +200,7 @@
void CohesiveTriaxialTest::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
+ body = shared_ptr<Body>(new Body); body->groupMask=2;
shared_ptr<CohFrictMat> physics(new CohFrictMat);
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
@@ -243,7 +243,7 @@
void CohesiveTriaxialTest::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents, bool wire)
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
+ body = shared_ptr<Body>(new Body); body->groupMask=2;
shared_ptr<CohFrictMat> physics(new CohFrictMat);
shared_ptr<Aabb> aabb(new Aabb);
=== modified file 'pkg/dem/FrictViscoPM.hpp'
--- pkg/dem/FrictViscoPM.hpp 2014-01-24 21:59:41 +0000
+++ pkg/dem/FrictViscoPM.hpp 2014-06-10 03:36:47 +0000
@@ -83,6 +83,7 @@
virtual void go(const shared_ptr<Material>& pp1, const shared_ptr<Material>& pp2, const shared_ptr<Interaction>& interaction);
FUNCTOR2D(FrictMat,FrictViscoMat);
+ DEFINE_FUNCTOR_ORDER_2D(FrictMat,FrictViscoMat);
YADE_CLASS_BASE_DOC_ATTRS(Ip2_FrictMat_FrictViscoMat_FrictViscoPhys,IPhysFunctor,"Converts a :yref:`FrictMat` and :yref:`FrictViscoMat` instance to :yref:`FrictViscoPhys` with corresponding parameters. Basically this functor corresponds to :yref:`Ip2_FrictMat_FrictMat_FrictPhys` with the only difference that damping in normal direction can be considered.",
((shared_ptr<MatchMaker>,kn,,,"Instance of :yref:`MatchMaker` determining how to compute interaction's normal contact stiffnesses. If this value is not given the elastic properties (i.e. young) of the two colliding materials are used to calculate the stiffness."))
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.cpp'
--- pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-05-28 08:47:32 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-06-10 14:27:43 +0000
@@ -270,11 +270,8 @@
//contactPhysics->kn = jointNormalStiffness*2.*R1*R2/(R1+R2); // very first expression from Luc
//contactPhysics->kn = (jkn1+jkn2)/2.0*2.*R1*R2/(R1+R2); // after putting jointNormalStiffness in material
- contactPhysics->kn = ( jkn1*Mathr::PI*pow(R1,2) + jkn2*Mathr::PI*pow(R2,2) )/2.0; // for a size independant expression
-
- //contactPhysics->ks = jointShearStiffness*2.*R1*R2/(R1+R2);
- //contactPhysics->ks = (jks1+jks2)/2.0*2.*R1*R2/(R1+R2);
- contactPhysics->ks = ( jks1*Mathr::PI*pow(R1,2) + jks2*Mathr::PI*pow(R2,2) )/2.0; // for a size independant expression
+ contactPhysics->kn = ( jkn1 + jkn2 ) /2.0 * contactPhysics->crossSection; // for a size independant expression
+ contactPhysics->ks = ( jks1 + jks2 ) /2.0 * contactPhysics->crossSection; // for a size independant expression
contactPhysics->tanDilationAngle = std::tan(std::min(jdil1,jdil2));
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.hpp'
--- pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-05-28 08:47:32 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-06-10 14:27:43 +0000
@@ -62,7 +62,7 @@
((bool,isOnJoint,false,,"defined as true when both interacting particles are :yref:`on joint<JCFpmState.onJoint>` and are in opposite sides of the joint surface. In this case, mechanical parameters of the interaction are derived from the ''joint...'' material properties of the particles, and normal of the interaction is re-oriented (see also :yref:`Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM`)."))
((Real,tanFrictionAngle,0,,"tangent of Coulomb friction angle for this interaction (auto. computed). [-]"))
((Real,crossSection,0,,"crossSection=pi*Rmin^2. [m2]"))
- ((Real,FnMax,0,,"computed from :yref:`tensile strength<JCFpmMat.tensileStrength>` (or joint variant) to define the maximum admissible normal force in traction. [N]"))
+ ((Real,FnMax,0,,"positiv value computed from :yref:`tensile strength<JCFpmMat.tensileStrength>` (or joint variant) to define the maximum admissible normal force in traction: Fn >= -FnMax. [N]"))
((Real,FsMax,0,,"computed from :yref:`cohesion<JCFpmMat.cohesion>` (or jointCohesion) to define the maximum admissible tangential force in shear, for Fn=0. [N]"))
((Vector3r,jointNormal,Vector3r::Zero(),,"normal direction to the joint, deduced from e.g. :yref:`<JCFpmState.jointNormal1>`."))
((Real,jointCumulativeSliding,0,,"sliding distance for particles interacting on a joint. Used, when :yref:`<Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM.smoothJoint>` is true, to take into account dilatancy due to shearing. [-]"))
=== modified file 'pkg/dem/SimpleShear.cpp'
--- pkg/dem/SimpleShear.cpp 2014-06-05 13:19:44 +0000
+++ pkg/dem/SimpleShear.cpp 2014-06-17 10:18:00 +0000
@@ -116,7 +116,7 @@
void SimpleShear::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius)
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(1);
+ body = shared_ptr<Body>(new Body); body->groupMask=1;
shared_ptr<NormalInelasticMat> mat(new NormalInelasticMat);
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
@@ -147,7 +147,7 @@
void SimpleShear::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents)
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(1);
+ body = shared_ptr<Body>(new Body); body->groupMask=1;
shared_ptr<NormalInelasticMat> mat(new NormalInelasticMat);
shared_ptr<Aabb> aabb(new Aabb);
// shared_ptr<BoxModel> gBox(new BoxModel);
=== modified file 'pkg/dem/SpheresFactory.cpp'
--- pkg/dem/SpheresFactory.cpp 2014-06-05 13:19:44 +0000
+++ pkg/dem/SpheresFactory.cpp 2014-06-17 10:18:00 +0000
@@ -182,7 +182,7 @@
b->shape=sphere;
b->state=state;
b->material=material;
- if (mask>0) {b->groupMask=Body::groupMask_t(mask);}
+ if (mask>0) {b->groupMask=mask;}
// insert particle in the simulation
scene->bodies->insert(b);
ids.push_back(b->getId());
=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp 2014-06-05 13:19:44 +0000
+++ pkg/dem/TesselationWrapper.cpp 2014-06-17 10:18:00 +0000
@@ -306,7 +306,7 @@
void createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool big, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
+ body = shared_ptr<Body>(new Body); body->groupMask=2;
shared_ptr<Sphere> iSphere(new Sphere);
body->state->blockedDOFs=State::DOF_NONE;
body->state->pos=position;
=== modified file 'pkg/dem/TriaxialTest.cpp'
--- pkg/dem/TriaxialTest.cpp 2014-06-05 13:19:44 +0000
+++ pkg/dem/TriaxialTest.cpp 2014-06-17 10:18:00 +0000
@@ -187,7 +187,7 @@
void TriaxialTest::createSphere(shared_ptr<Body>& body, Vector3r position, Real radius, bool big, bool dynamic )
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
+ body = shared_ptr<Body>(new Body); body->groupMask=2;
shared_ptr<Aabb> aabb(new Aabb);
shared_ptr<Sphere> iSphere(new Sphere);
body->state->blockedDOFs=State::DOF_NONE;
@@ -213,7 +213,7 @@
void TriaxialTest::createBox(shared_ptr<Body>& body, Vector3r position, Vector3r extents, bool wire)
{
- body = shared_ptr<Body>(new Body); body->groupMask=Body::groupMask_t(2);
+ body = shared_ptr<Body>(new Body); body->groupMask=2;
body->state->blockedDOFs=State::DOF_ALL;
shared_ptr<Aabb> aabb(new Aabb);
aabb->color = Vector3r(1,0,0);
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2014-06-05 13:19:44 +0000
+++ pkg/dem/VTKRecorder.cpp 2014-06-17 10:18:00 +0000
@@ -37,8 +37,8 @@
YADE_PLUGIN((VTKRecorder));
CREATE_LOGGER(VTKRecorder);
-#ifdef BODY_GROUP_MASK_ARBITRARY_PRECISION
-#define GET_MASK(b) boost::python::extract<int>(b->groupMask)
+#ifdef YADE_MASK_ARBITRARY
+#define GET_MASK(b) b->groupMask.to_ulong()
#else
#define GET_MASK(b) b->groupMask
#endif
@@ -479,10 +479,12 @@
spheresCoordNumbSPH->InsertNextValue(b->coordNumber());
#endif
#ifdef YADE_LIQMIGRATION
- spheresLiqVol->InsertNextValue(b->Vf);
- const Real tmpVolIter = liqVolIterBody(b);
- spheresLiqVolIter->InsertNextValue(tmpVolIter);
- spheresLiqVolTotal->InsertNextValue(tmpVolIter + b->Vf);
+ if (recActive[REC_LIQ]) {
+ spheresLiqVol->InsertNextValue(b->Vf);
+ const Real tmpVolIter = liqVolIterBody(b);
+ spheresLiqVolIter->InsertNextValue(tmpVolIter);
+ spheresLiqVolTotal->InsertNextValue(tmpVolIter + b->Vf);
+ }
#endif
if (recActive[REC_MATERIALID]) spheresMaterialId->InsertNextValue(b->material->id);
continue;
@@ -626,9 +628,11 @@
spheresUg->GetPointData()->AddArray(spheresCoordNumbSPH);
#endif
#ifdef YADE_LIQMIGRATION
- spheresUg->GetPointData()->AddArray(spheresLiqVol);
- spheresUg->GetPointData()->AddArray(spheresLiqVolIter);
- spheresUg->GetPointData()->AddArray(spheresLiqVolTotal);
+ if (recActive[REC_LIQ]) {
+ spheresUg->GetPointData()->AddArray(spheresLiqVol);
+ spheresUg->GetPointData()->AddArray(spheresLiqVolIter);
+ spheresUg->GetPointData()->AddArray(spheresLiqVolTotal);
+ }
#endif
if (recActive[REC_STRESS]){
spheresUg->GetPointData()->AddArray(spheresNormalStressVec);
=== modified file 'pkg/dem/VTKRecorder.hpp'
--- pkg/dem/VTKRecorder.hpp 2014-05-08 05:57:04 +0000
+++ pkg/dem/VTKRecorder.hpp 2014-06-12 15:11:52 +0000
@@ -17,7 +17,7 @@
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.",
((bool,compress,false,,"Compress output XML files [experimental]."))
((bool,ascii,false,,"Store data as readable text in the XML file (sets `vtkXMLWriter <http://www.vtk.org/doc/nightly/html/classvtkXMLWriter.html>`__ data mode to ``vtkXMLWriter::Ascii``, while the default is ``Appended``"))
- ((bool,skipFacetIntr,true,,"Skip interactions with facets, when saving interactions"))
+ ((bool,skipFacetIntr,true,,"Skip interactions that are not of sphere-sphere type (e.g. sphere-facet, sphere-box...), when saving interactions"))
((bool,skipNondynamic,false,,"Skip non-dynamic spheres (but not facets)."))
#ifdef YADE_VTK_MULTIBLOCK
((bool,multiblock,false,,"Use multi-block (``.vtm``) files to store data, rather than separate ``.vtu`` files."))
=== modified file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 2014-05-27 11:24:16 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-06-18 07:35:57 +0000
@@ -361,8 +361,13 @@
// Calculate, how much new contacts will be at each body
for (unsigned int i=0; i<scene->addIntrs.size(); i++) {
- addBodyMapInt( bI, scene->addIntrs[i]->getId1() );
- addBodyMapInt( bI, scene->addIntrs[i]->getId2() );
+ shared_ptr<Body> b1 = Body::byId(scene->addIntrs[i]->getId1(),scene);
+ shared_ptr<Body> b2 = Body::byId(scene->addIntrs[i]->getId2(),scene);
+
+ if(not(mask!=0 && ((b1->groupMask & b2->groupMask & mask)==0))) {
+ addBodyMapInt( bI, scene->addIntrs[i]->getId1() );
+ addBodyMapInt( bI, scene->addIntrs[i]->getId2() );
+ }
}
// Update volume water at each deleted interaction for each body
@@ -384,43 +389,42 @@
const id_t id2 = b2->id;
ViscElCapPhys* Vb=dynamic_cast<ViscElCapPhys*>(scene->addIntrs[i]->phys.get());
- const Real Vmax = vMax(b1, b2);
- Vb->Vmax = Vmax;
-
- Real Vf1 = 0.0;
- Real Vf2 = 0.0;
-
- if ((b1->Vmin)<b1->Vf) {
- Vf1 = (b1->Vf - b1->Vmin)/bI[id1];
- }
-
- if ((b2->Vmin)<b2->Vf) {
- Vf2 = (b2->Vf - b2->Vmin)/bI[id2];
- }
-
- Real Vrup = Vf1+Vf2;
-
+
if(mask!=0 && ((b1->groupMask & b2->groupMask & mask)==0)) {
- Vf1 = 0;
- Vf2 = 0;
- Vrup = 0;
- } else if (Vrup > Vmax) {
- Vf1 *= Vmax/Vrup;
- Vf2 *= Vmax/Vrup;
- Vrup = Vf1 + Vf2;
- }
-
- liqVolShr += Vrup;
- addBodyMapReal(bodyUpdateLiquid, id1, -Vf1);
- addBodyMapReal(bodyUpdateLiquid, id2, -Vf2);
-
- Vb->Vb = Vrup;
- if (particleconserve) {
- Vb->Vf1 = Vf1;
- Vb->Vf2 = Vf2;
+ Vb->Vb = 0.0;
+ Vb->Vf1 = 0.0;
+ Vb->Vf2 = 0.0;
+ Vb->Vmax = 0.0;
} else {
- Vb->Vf1 = Vrup/2.0;
- Vb->Vf2 = Vrup/2.0;
+ const Real Vmax = vMax(b1, b2);
+ Vb->Vmax = Vmax;
+
+ Real Vf1 = 0.0;
+ Real Vf2 = 0.0;
+
+ if ((b1->Vmin)<b1->Vf) { Vf1 = (b1->Vf - b1->Vmin)/bI[id1]; }
+ if ((b2->Vmin)<b2->Vf) { Vf2 = (b2->Vf - b2->Vmin)/bI[id2]; }
+
+ Real Vrup = Vf1+Vf2;
+
+ if (Vrup > Vmax) {
+ Vf1 *= Vmax/Vrup;
+ Vf2 *= Vmax/Vrup;
+ Vrup = Vf1 + Vf2;
+ }
+
+ liqVolShr += Vrup;
+ addBodyMapReal(bodyUpdateLiquid, id1, -Vf1);
+ 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;
+ }
}
}
@@ -435,7 +439,6 @@
for (mapBodyInt::const_iterator it = bodyNeedUpdate.begin(); it != bodyNeedUpdate.end(); ++it) {
updateLiquid(Body::byId(it->first));
}
-
}
void LiqControl::updateLiquid(shared_ptr<Body> b){
@@ -451,7 +454,7 @@
for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
if(!((*it).second) or !(((*it).second)->isReal())) continue;
ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
- if (physT->Vb<physT->Vmax) {
+ if ((physT->Vb < physT->Vmax) and (physT->Vmax > 0)) {
LiqContactsAccept+=physT->Vmax-physT->Vb;
contactN++;
}
@@ -472,15 +475,15 @@
for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
if(!((*it).second) or !(((*it).second)->isReal())) continue;
ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
- if (physT->Vb<physT->Vmax) {
+ if ((physT->Vb < physT->Vmax) and (physT->Vmax > 0)) {
const Real addVolLiq = (physT->Vmax - physT->Vb)*FillLevel;
liqVolShr += addVolLiq;
physT->Vb += addVolLiq;
if (particleconserve) {
if (((*it).second)->getId1() == (*it).first) {
+ physT->Vf2+=addVolLiq;
+ } else if (((*it).second)->getId2() == (*it).first) {
physT->Vf1+=addVolLiq;
- } else if (((*it).second)->getId2() == (*it).first) {
- physT->Vf2+=addVolLiq;
}
} else {
physT->Vf1+=addVolLiq/2.0;
@@ -531,14 +534,30 @@
Real liqVolIterBody (shared_ptr<Body> b) {
Real LiqVol = 0.0;
- for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
- if(!((*it).second) or !(((*it).second)->isReal())) continue;
- ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
- if (physT->Vb>0) {
- LiqVol += physT->Vb/2.0;
+ if (!b) {
+ return LiqVol;
+ } else {
+ for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
+ if(!((*it).second) or !(((*it).second)->isReal())) continue;
+ ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
+ if (physT and physT->Vb and physT->Vb>0) {
+ if (((*it).second)->id1 == b->id) {
+ if (physT->Vf1 > 0 or physT->Vf2 > 0) {
+ LiqVol += physT->Vf1;
+ } else {
+ LiqVol += physT->Vb/2.0;
+ }
+ } else {
+ if (physT->Vf1 > 0 or physT->Vf2 > 0) {
+ LiqVol += physT->Vf2;
+ } else {
+ LiqVol += physT->Vb/2.0;
+ }
+ }
+ }
}
+ return LiqVol;
}
- return LiqVol;
}
Real LiqControl::liqVolBody (id_t id) const {
@@ -564,4 +583,22 @@
}
return totalLiqVol;
}
+
+bool LiqControl::addLiqInter(id_t id1, id_t id2, Real liq) {
+ if (id1<0 or id2<0 or id1==id2 or liq<=0) return false;
+
+ Scene* scene=Omega::instance().getScene().get();
+ shared_ptr<InteractionContainer>& intrs=scene->interactions;
+ const shared_ptr<Interaction>& I=intrs->find(id1,id2);
+ if (I->isReal()) {
+ ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(I->phys.get());
+ if (physT and physT->Vb <= physT->Vmax and liq <= (physT->Vmax - physT->Vb)) {
+ physT->Vb += liq;
+ physT->Vf1 += liq/2.0;
+ physT->Vf2 += liq/2.0;
+ return true;
+ }
+ }
+ return false;
+}
#endif
=== modified file 'pkg/dem/ViscoelasticCapillarPM.hpp'
--- pkg/dem/ViscoelasticCapillarPM.hpp 2014-05-27 14:20:39 +0000
+++ pkg/dem/ViscoelasticCapillarPM.hpp 2014-06-13 13:58:12 +0000
@@ -92,6 +92,7 @@
Real vMax(shared_ptr<Body> b1, shared_ptr<Body> b2);
Real totalLiqVol(int mask) const;
Real liqVolBody(id_t id) const;
+ bool addLiqInter(id_t id1, id_t id2, Real liq);
void updateLiquid(shared_ptr<Body> b);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(LiqControl,PartialEngine,"This engine implements liquid migration model, introduced here [Mani2013]_ . ",
((int,mask,0,, "Bitmask for liquid creation."))
@@ -103,6 +104,7 @@
,/* py */
.def("totalLiq",&LiqControl::totalLiqVol,(boost::python::arg("mask")=0),"Return total volume of water in simulation.")
.def("liqBody",&LiqControl::liqVolBody,(boost::python::arg("id")=-1),"Return total volume of water in body.")
+ .def("addLiqInter",&LiqControl::addLiqInter,(boost::python::arg("id1")=-1, boost::python::arg("id2")=-1, boost::python::arg("liq")=-1),"Add liquid into the interaction.")
);
};
=== modified file 'pkg/pfv/DFNFlow.cpp'
--- pkg/pfv/DFNFlow.cpp 2014-05-27 10:33:52 +0000
+++ pkg/pfv/DFNFlow.cpp 2014-06-20 20:10:04 +0000
@@ -11,11 +11,11 @@
//it will save compilation time for everyone else
//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
// #define DFNFLOW
+
#ifdef DFNFLOW
-#define TEMPLATE_FLOW_NAME DFNFlowEngineT
-#include <yade/pkg/pfv/FlowEngine.hpp>
+#include "FlowEngine_DFNFlowEngineT.hpp"
-class DFNCellInfo : public FlowCellInfo
+class DFNCellInfo : public FlowCellInfo_DFNFlowEngineT
{
public:
Real anotherVariable;
@@ -23,12 +23,12 @@
};
-class DFNVertexInfo : public FlowVertexInfo {
+class DFNVertexInfo : public FlowVertexInfo_DFNFlowEngineT {
public:
//same here if needed
};
-typedef TemplateFlowEngine<DFNCellInfo,DFNVertexInfo> DFNFlowEngineT;
+typedef TemplateFlowEngine_DFNFlowEngineT<DFNCellInfo,DFNVertexInfo> DFNFlowEngineT;
REGISTER_SERIALIZABLE(DFNFlowEngineT);
YADE_PLUGIN((DFNFlowEngineT));
@@ -104,4 +104,4 @@
#endif //DFNFLOW
-#endif //FLOWENGINE
\ No newline at end of file
+#endif //FLOWENGINE
=== modified file 'pkg/pfv/DummyFlowEngine.cpp'
--- pkg/pfv/DummyFlowEngine.cpp 2014-05-23 13:20:43 +0000
+++ pkg/pfv/DummyFlowEngine.cpp 2014-06-20 15:28:21 +0000
@@ -13,27 +13,27 @@
//when you want it compiled, you can pass -DDFNFLOW to cmake, or just uncomment the following line
// #define DUMMYFLOW
#ifdef DUMMYFLOW
-#define TEMPLATE_FLOW_NAME DummyFlowEngineT
-#include <yade/pkg/pfv/FlowEngine.hpp>
+
+#include "FlowEngine_DummyFlowEngineT.hpp"
/// We can add data to the Info types by inheritance
-class DummyCellInfo : public FlowCellInfo
+class DummyCellInfo : public FlowCellInfo_DummyFlowEngineT
{
public:
Real anotherVariable;
void anotherFunction() {};
};
-class DummyVertexInfo : public FlowVertexInfo {
+class DummyVertexInfo : public FlowVertexInfo_DummyFlowEngineT {
public:
//same here if needed
};
-typedef TemplateFlowEngine<DummyCellInfo,DummyVertexInfo> TEMPLATE_FLOW_NAME;
-REGISTER_SERIALIZABLE(TEMPLATE_FLOW_NAME);
-YADE_PLUGIN((TEMPLATE_FLOW_NAME));
+typedef TemplateFlowEngine_DummyFlowEngineT<DummyCellInfo,DummyVertexInfo> DummyFlowEngineT;
+REGISTER_SERIALIZABLE(DummyFlowEngineT);
+YADE_PLUGIN((DummyFlowEngineT));
-class DummyFlowEngine : public TEMPLATE_FLOW_NAME
+class DummyFlowEngine : public DummyFlowEngineT
{
public :
//We can overload every functions of the base engine to make it behave differently
@@ -42,9 +42,9 @@
//If a new function is specific to the derived engine, put it here, else go to the base TemplateFlowEngine
//if it is useful for everyone
- void fancyFunction(Real what); {cerr<<"yes, I'm a new function"<<end;}
+ void fancyFunction(Real what);
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DummyFlowEngine,TEMPLATE_FLOW_NAME,"documentation here",
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(DummyFlowEngine,DummyFlowEngineT,"documentation here",
((Real, myNewAttribute, 0,,"useless example"))
,/*DummyFlowEngineT()*/,
,
@@ -55,6 +55,5 @@
REGISTER_SERIALIZABLE(DummyFlowEngine);
YADE_PLUGIN((DummyFlowEngine));
-void DummyFlowEngine::fancyFunction(Real what) {cerr<<"yes, I'm a new function"<<end;}
-#undef TEMPLATE_FLOW_NAME DummyFlowEngineT //To be sure it will not conflict, maybe not needed
+void DummyFlowEngine::fancyFunction(Real what) {std::cerr<<"yes, I'm a new function"<<std::endl;}
#endif //DummyFLOW
=== modified file 'pkg/pfv/FlowEngine.cpp'
--- pkg/pfv/FlowEngine.cpp 2014-04-03 13:37:55 +0000
+++ pkg/pfv/FlowEngine.cpp 2014-06-20 17:48:57 +0000
@@ -8,13 +8,12 @@
#ifdef YADE_CGAL
#ifdef FLOW_ENGINE
-#define TEMPLATE_FLOW_NAME FlowEngineT
-#include "FlowEngine.hpp"
+#include "FlowEngine_FlowEngineT.hpp"
// To register properly, we need to first instantiate an intermediate class, then inherit from it with correct class names in YADE_CLASS macro
// The intermediate one would be seen with the name "TemplateFlowEngine" by python, thus it would not work when more than one class are derived, they would all
// be named "TemplateFlowEngine" ...
-typedef TemplateFlowEngine<FlowCellInfo,FlowVertexInfo> FlowEngineT;
+typedef TemplateFlowEngine_FlowEngineT<FlowCellInfo_FlowEngineT,FlowVertexInfo_FlowEngineT> FlowEngineT;
REGISTER_SERIALIZABLE(FlowEngineT);
class FlowEngine : public FlowEngineT
=== modified file 'pkg/pfv/FlowEngine.hpp'
--- pkg/pfv/FlowEngine.hpp 2014-06-04 17:19:50 +0000
+++ pkg/pfv/FlowEngine.hpp 2014-06-20 15:28:21 +0000
@@ -1,45 +1,4 @@
-/*************************************************************************
-* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
-* Copyright (C) 2009 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. *
-*************************************************************************/
-
-/*!
-
-FlowEngine is an interface between Yade and the fluid solvers defined in lib/triangulation and using the PFV scheme.
-There are also a few non-trivial functions defined here, such has building triangulation and computating elements volumes.
-The code strongly relies on CGAL library for building triangulations on the top of sphere packings.
-CGAL's trangulations introduce the data associated to each element of the triangulation through template parameters (Cell_info and Vertex_info).
-FlowEngine is thus defined via the template class TemplateFlowEngine(Cellinfo,VertexInfo), for easier addition of variables in various couplings.
-PeriodicFlowEngine is a variant for periodic boundary conditions.
-
-Which solver will be actually used internally to obtain pore pressure will depend partly on compile time flags (libcholmod available and #define LINSOLV),
-and on runtime settings (useSolver=0: Gauss-Seidel (iterative), useSolver=1: CHOLMOD if available (direct sparse cholesky)).
-
-The files defining lower level classes are in yade/lib. The code uses CGAL::Triangulation3 for managing the mesh and storing data. Eigen3::Sparse, suitesparse::cholmod, and metis are used for solving the linear systems with a direct method (Cholesky). An iterative method (Gauss-Seidel) is implemented directly in Yade and can be used as a fallback (see FlowEngine::useSolver).
-
-Most classes in lib/triangulation are templates, and are therefore completely defined in header files.
-A pseudo hpp/cpp split is reproduced for clarity, with hpp/ipp extensions (but again, in the end they are all include files).
-The files are
-- RegularTriangulation.h : declaration of info types embeded in the mesh (template parameters of CGAL::Triangulation3), and instanciation of RTriangulation classes
-- Tesselation.h/ipp : encapsulate RegularTriangulation and adds functions to manipulate the dual (Voronoi) graph of the triangulation
-- Network.hpp/ipp : specialized for PFV model (the two former are used independently by TesselationWrapper), a set of functions to determine volumes and surfaces of intersections between spheres and various subdomains. Contains two triangulations for smooth transitions while remeshing - e.g. interpolating values in the new mesh using the previous one.
-- FlowBoundingSphere.hpp/ipp and PeriodicFlow.hpp/ipp + LinSolv variants: implement the solver in itself (mesh, boundary conditions, solving, defining fluid-particles interactions)
-- FlowEngine.hpp/ipp/cpp (this file)
-
-Variants for periodic boundary conditions are also present.
-
-*/
-
#pragma once
-#include<yade/core/PartialEngine.hpp>
-#include<yade/pkg/dem/TriaxialCompressionEngine.hpp>
-#include<yade/pkg/dem/TesselationWrapper.hpp>
-#include<yade/lib/triangulation/FlowBoundingSphere.hpp>
-#include<yade/lib/triangulation/PeriodicFlow.hpp>
-
/// Frequently used:
typedef CGT::CVector CVector;
typedef CGT::Point Point;
@@ -50,391 +9,8 @@
inline Vector3r makeVector3r ( const Point& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
inline Vector3r makeVector3r ( const CVector& yv ) {return Vector3r ( yv[0],yv[1],yv[2] );}
-/// C++ templates and YADE_CLASS_... macro family are not very compatible, this #define is a hack to make it work
-/// TEMPLATE_FLOW_NAME will be the name of a serializable TemplateFlowEngine<...> instance, which can in turn be
-/// inherited from. The instance itself will be useless for actual simulations.
-#ifndef TEMPLATE_FLOW_NAME
-#error You must define TEMPLATE_FLOW_NAME in your source file before including FlowEngine.hpp
-#endif
-
#ifdef LINSOLV
#define DEFAULTSOLVER CGT::FlowBoundingSphereLinSolv<_Tesselation>
#else
#define DEFAULTSOLVER CGT::FlowBoundingSphere<_Tesselation>
#endif
-
-template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=DEFAULTSOLVER >
-class TemplateFlowEngine : public PartialEngine
-{
- public :
- typedef solverT FlowSolver;
- typedef FlowSolver Solver;//FIXME: useless alias, search/replace then remove
- typedef typename FlowSolver::Tesselation Tesselation;
- typedef typename FlowSolver::RTriangulation RTriangulation;
- typedef typename FlowSolver::FiniteVerticesIterator FiniteVerticesIterator;
- typedef typename FlowSolver::FiniteCellsIterator FiniteCellsIterator;
- typedef typename FlowSolver::CellHandle CellHandle;
- typedef typename FlowSolver::FiniteEdgesIterator FiniteEdgesIterator;
- typedef typename FlowSolver::VertexHandle VertexHandle;
- typedef typename RTriangulation::Triangulation_data_structure::Cell::Info CellInfo;
- typedef typename RTriangulation::Triangulation_data_structure::Vertex::Info VertexInfo;
-
-// protected:
- shared_ptr<FlowSolver> solver;
- shared_ptr<FlowSolver> backgroundSolver;
- volatile bool backgroundCompleted;
- Cell cachedCell;
- struct posData {Body::id_t id; Vector3r pos; Real radius; bool isSphere; bool exists; posData(){exists=0;}};
- 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
- virtual void setPositionsBuffer(bool current);
- virtual void trickPermeability() {};
-
- public :
- int retriangulationLastIter;
- enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax};
- Vector3r normal [6];
- bool currentTes;
- bool metisForced;
- int idOffset;
- double epsVolCumulative;
- int ReTrg;
- int ellapsedIter;
- void initSolver (FlowSolver& flow);
- #ifdef LINSOLV
- void setForceMetis (bool force);
- bool getForceMetis ();
- #endif
- void triangulate (Solver& flow);
- void addBoundary (Solver& flow);
- void buildTriangulation (double pZero, Solver& flow);
- void buildTriangulation (Solver& flow);
- void updateVolumes (Solver& flow);
- void initializeVolumes (Solver& flow);
- void boundaryConditions(Solver& flow);
- void updateBCs () {
- if (solver->tesselation().maxId>0) boundaryConditions(*solver);//avoids crash at iteration 0, when the packing is not bounded yet
- else LOG_ERROR("updateBCs not applied");
- solver->pressureChanged=true;}
-
- void imposeFlux(Vector3r pos, Real flux);
- unsigned int imposePressure(Vector3r pos, Real p);
- void setImposedPressure(unsigned int cond, Real p);
- void clearImposedPressure();
- void clearImposedFlux();
- void computeViscousForces(Solver& flow);
- Real getCellFlux(unsigned int cond);
- Real getBoundaryFlux(unsigned int boundary) {return solver->boundaryFlux(boundary);}
- Vector3r fluidForce(unsigned int idSph) {
- const CGT::CVector& f=solver->T[solver->currentTes].vertex(idSph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
- Vector3r averageVelocity();
-
- Vector3r shearLubForce(unsigned int id_sph) {
- return (solver->shearLubricationForces.size()>id_sph)?solver->shearLubricationForces[id_sph]:Vector3r::Zero();}
- Vector3r shearLubTorque(unsigned int id_sph) {
- return (solver->shearLubricationTorques.size()>id_sph)?solver->shearLubricationTorques[id_sph]:Vector3r::Zero();}
- Vector3r pumpLubTorque(unsigned int id_sph) {
- return (solver->pumpLubricationTorques.size()>id_sph)?solver->pumpLubricationTorques[id_sph]:Vector3r::Zero();}
- Vector3r twistLubTorque(unsigned int id_sph) {
- return (solver->twistLubricationTorques.size()>id_sph)?solver->twistLubricationTorques[id_sph]:Vector3r::Zero();}
- Vector3r normalLubForce(unsigned int id_sph) {
- return (solver->normalLubricationForce.size()>id_sph)?solver->normalLubricationForce[id_sph]:Vector3r::Zero();}
- Matrix3r bodyShearLubStress(unsigned int id_sph) {
- return (solver->shearLubricationBodyStress.size()>id_sph)?solver->shearLubricationBodyStress[id_sph]:Matrix3r::Zero();}
- Matrix3r bodyNormalLubStress(unsigned int id_sph) {
- return (solver->normalLubricationBodyStress.size()>id_sph)?solver->normalLubricationBodyStress[id_sph]:Matrix3r::Zero();}
- Vector3r shearVelocity(unsigned int interaction) {
- return (solver->deltaShearVel[interaction]);}
- Vector3r normalVelocity(unsigned int interaction) {
- return (solver->deltaNormVel[interaction]);}
- Matrix3r normalStressInteraction(unsigned int interaction) {
- return (solver->normalStressInteraction[interaction]);}
- Matrix3r shearStressInteraction(unsigned int interaction) {
- return (solver->shearStressInteraction[interaction]);}
- Vector3r normalVect(unsigned int interaction) {
- return (solver->normalV[interaction]);}
- Real surfaceDistanceParticle(unsigned int interaction) {
- return (solver->surfaceDistance[interaction]);}
- Real edgeSize() {
- return (solver->edgeIds.size());}
- Real OSI() {
- return (solver->onlySpheresInteractions.size());}
- int onlySpheresInteractions(unsigned int interaction) {
- return (solver->onlySpheresInteractions[interaction]);}
- boost::python::list getConstrictions(bool all) {
- vector<Real> csd=solver->getConstrictions(); boost::python::list pycsd;
- for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k]!=0) || csd[k]>0) pycsd.append(csd[k]); return pycsd;}
- boost::python::list getConstrictionsFull(bool all) {
- vector<Constriction> csd=solver->getConstrictionsFull(); boost::python::list pycsd;
- for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k].second[0]!=0) || csd[k].second[0]>0) {
- boost::python::list cons;
- cons.append(csd[k].first.first);
- cons.append(csd[k].first.second);
- cons.append(csd[k].second[0]);
- cons.append(csd[k].second[1]);
- cons.append(csd[k].second[2]);
- cons.append(csd[k].second[3]);
- pycsd.append(cons);}
- return pycsd;}
-
- template<class Cellhandle>
- Real volumeCellSingleFictious (Cellhandle cell);
- template<class Cellhandle>
- Real volumeCellDoubleFictious (Cellhandle cell);
- template<class Cellhandle>
- Real volumeCellTripleFictious (Cellhandle cell);
- template<class Cellhandle>
- Real volumeCell (Cellhandle cell);
- void Oedometer_Boundary_Conditions();
- void averageRealCellVelocity();
- void saveVtk(const char* folder) {solver->saveVtk(folder);}
- vector<Real> avFlVelOnSph(unsigned int idSph) {return solver->averageFluidVelocityOnSphere(idSph);}
-
-// void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; updateTriangulation=true;}
- void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}
- double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
- int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}
- unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}
- boost::python::list getVertices(unsigned int id){
- boost::python::list ids;
- if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
- for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->vertex(i)->info().id());
- return ids;
- }
- double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
- double averagePressure(){return solver->averagePressure();}
-
- void emulateAction(){
- scene = Omega::instance().getScene().get();
- action();}
- #ifdef LINSOLV
- void exportMatrix(string filename) {if (useSolver==3) solver->exportMatrix(filename.c_str());
- else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
- void exportTriplets(string filename) {if (useSolver==3) solver->exportTriplets(filename.c_str());
- else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
- void cholmodStats() {
- cerr << cholmod_print_common((char*)string("PFV Cholmod factorization").c_str(),&(solver->eSolver.cholmod()))<<endl;
- cerr << "cholmod method:" << solver->eSolver.cholmod().selected<<endl;
- cerr << "METIS called:"<<solver->eSolver.cholmod().called_nd<<endl;}
- bool metisUsed() {return bool(solver->eSolver.cholmod().called_nd);}
- #endif
-
- virtual ~TemplateFlowEngine();
- virtual void action();
- virtual void backgroundAction();
-
- //commodities
- void compTessVolumes() {
- solver->T[solver->currentTes].compute();
- solver->T[solver->currentTes].computeVolumes();
- }
- Real getVolume (Body::id_t id) {
- if (solver->T[solver->currentTes].Max_id() <= 0) {emulateAction(); /*LOG_WARN("Not triangulated yet, emulating action");*/}
- if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes();/* LOG_WARN("Computing all volumes now, as you did not request it explicitely.");*/}
- return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
-
- YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine,TEMPLATE_FLOW_NAME,PartialEngine,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. See the example script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n.. note::Multi-threading seems to work fine for Cholesky decomposition, but it fails for the solve phase in which -j1 is the fastest, here we specify thread numbers independently using :yref:`FlowEngine::numFactorizeThreads` and :yref:`FlowEngine::numSolveThreads`. These multhreading settings are only impacting the behaviour of openblas library and are relatively independant of :yref:`FlowEngine::multithread`. However, the settings have to be globally consistent. For instance, :yref:`multithread<FlowEngine::multithread>` =True with yref:`numFactorizeThreads<FlowEngine::numFactorizeThreads>` = yref:`numSolveThreads<FlowEngine::numSolveThreads>` = 4 implies that openblas will mobilize 8 processors at some point. If the system does not have so many procs. it will hurt performance.",
- ((bool,isActivated,true,,"Activates Flow Engine"))
- ((bool,first,true,,"Controls the initialization/update phases"))
- ((bool,doInterpolate,false,,"Force the interpolation of cell's info while remeshing. By default, interpolation would be done only for compressible fluids. It can be forced with this flag."))
- ((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
- ((Real, dt, 0,,"timestep [s]"))
- ((bool,permeabilityMap,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
- ((bool, slipBoundary, true,, "Controls friction condition on lateral walls"))
- ((bool,waveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
- ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))
- ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
- ((bool, debug, false,,"Activate debug messages"))
- ((double, wallThickness,0.001,,"Walls thickness"))
- ((double,pZero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
- ((double,tolerance,1e-06,,"Gauss-Seidel tolerance"))
- ((double,relax,1.9,,"Gauss-Seidel relaxation"))
- ((bool, updateTriangulation, 0,,"If true the medium is retriangulated. Can be switched on to force retriangulation after some events (else it will be true periodicaly based on :yref:`FlowEngine::defTolerance` and :yref:`FlowEngine::meshUpdateInterval`. Of course, it costs CPU time."))
- ((int,meshUpdateInterval,1000,,"Maximum number of timesteps between re-triangulation events. See also :yref:`FlowEngine::defTolerance`."))
- ((double, epsVolMax, 0,(Attr::readonly),"Maximal absolute volumetric strain computed at each iteration. |yupdate|"))
- ((double, defTolerance,0.05,,"Cumulated deformation threshold for which retriangulation of pore space is performed. If negative, the triangulation update will occure with a fixed frequency on the basis of :yref:`FlowEngine::meshUpdateInterval`"))
- ((double, porosity, 0,(Attr::readonly),"Porosity computed at each retriangulation |yupdate|"))
- ((bool,meanKStat,false,,"report the local permeabilities' correction"))
- ((bool,clampKValues,true,,"If true, clamp local permeabilities in [minKdivKmean,maxKdivKmean]*globalK. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))
- ((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
- ((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
- ((double,permeabilityFactor,1.0,,"permability multiplier"))
- ((double,viscosity,1.0,,"viscosity of the fluid"))
- ((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))
- ((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
- ((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))
- ((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, ymax,3,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, zmin,4,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
- ((int, zmax,5,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
-
- ((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends."))
- ((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future)."))
- //FIXME: to be implemented:
- ((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
- ((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
- ((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))
- ((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
-// ((bool, display_force, false,,"display the lubrication force applied on particles"))
-// ((bool, create_file, false,,"create file of velocities"))
- ((bool, viscousShear, false,,"compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))
- ((bool, shearLubrication, false,,"compute shear lubrication force as developped by Brule (FIXME: ref.) "))
- ((bool, pumpTorque, false,,"Compute pump torque applied on particles "))
- ((bool, twistTorque, false,,"Compute twist torque applied on particles "))
- ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
- ((bool, pressureForce, true,,"compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
- ((bool, normalLubrication, false,,"compute normal lubrication force as developped by Brule"))
- ((bool, viscousNormalBodyStress, false,,"compute normal viscous stress applied on each body"))
- ((bool, viscousShearBodyStress, false,,"compute shear viscous stress applied on each body"))
- ((bool, multithread, false,,"Build triangulation and factorize in the background (multi-thread mode)"))
- #ifdef EIGENSPARSE_LIB
- ((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
- ((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
- #endif
- ((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`TEMPLATE_FLOW_NAME::boundaryXPos`"))
- ((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`TEMPLATE_FLOW_NAME::boundaryPressure`"))
- ,
- /*deprec*/
- ((meanK_opt,clampKValues,"the name changed"))
- ,,
- timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
- for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero(); wallIds[i]=i;}
- normal[wall_ymin].y()=normal[wall_xmin].x()=normal[wall_zmin].z()=1;
- normal[wall_ymax].y()=normal[wall_xmax].x()=normal[wall_zmax].z()=-1;
- solver = shared_ptr<FlowSolver> (new FlowSolver);
- first=true;
- epsVolMax=epsVolCumulative=retriangulationLastIter=0;
- ReTrg=1;
- backgroundCompleted=true;
- ellapsedIter=0;
- metisForced=false;
- ,
- .def("imposeFlux",&TemplateFlowEngine::imposeFlux,(boost::python::arg("pos"),boost::python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
- .def("imposePressure",&TemplateFlowEngine::imposePressure,(boost::python::arg("pos"),boost::python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
- .def("setImposedPressure",&TemplateFlowEngine::setImposedPressure,(boost::python::arg("cond"),boost::python::arg("p")),"Set pressure value at the point indexed 'cond'.")
- .def("clearImposedPressure",&TemplateFlowEngine::clearImposedPressure,"Clear the list of points with pressure imposed.")
- .def("clearImposedFlux",&TemplateFlowEngine::clearImposedFlux,"Clear the list of points with flux imposed.")
- .def("getCellFlux",&TemplateFlowEngine::getCellFlux,(boost::python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
- .def("getBoundaryFlux",&TemplateFlowEngine::getBoundaryFlux,(boost::python::arg("boundary")),"Get total flux through boundary defined by its body id.\n\n.. note:: The flux may be not zero even for no-flow condition. This artifact comes from cells which are incident to two or more boundaries (along the edges of the sample, typically). Such flux evaluation on impermeable boundary is just irrelevant, it does not imply that the boundary condition is not applied properly.")
- .def("getConstrictions",&TemplateFlowEngine::getConstrictions,(boost::python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
- .def("getConstrictionsFull",&TemplateFlowEngine::getConstrictionsFull,(boost::python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
- .def("edgeSize",&TemplateFlowEngine::edgeSize,"Return the number of interactions.")
- .def("OSI",&TemplateFlowEngine::OSI,"Return the number of interactions only between spheres.")
-
- .def("saveVtk",&TemplateFlowEngine::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
- .def("avFlVelOnSph",&TemplateFlowEngine::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")
- .def("fluidForce",&TemplateFlowEngine::fluidForce,(boost::python::arg("idSph")),"Return the fluid force on sphere idSph.")
- .def("shearLubForce",&TemplateFlowEngine::shearLubForce,(boost::python::arg("idSph")),"Return the shear lubrication force on sphere idSph.")
- .def("shearLubTorque",&TemplateFlowEngine::shearLubTorque,(boost::python::arg("idSph")),"Return the shear lubrication torque on sphere idSph.")
- .def("normalLubForce",&TemplateFlowEngine::normalLubForce,(boost::python::arg("idSph")),"Return the normal lubrication force on sphere idSph.")
- .def("bodyShearLubStress",&TemplateFlowEngine::bodyShearLubStress,(boost::python::arg("idSph")),"Return the shear lubrication stress on sphere idSph.")
- .def("bodyNormalLubStress",&TemplateFlowEngine::bodyNormalLubStress,(boost::python::arg("idSph")),"Return the normal lubrication stress on sphere idSph.")
- .def("shearVelocity",&TemplateFlowEngine::shearVelocity,(boost::python::arg("idSph")),"Return the shear velocity of the interaction.")
- .def("normalVelocity",&TemplateFlowEngine::normalVelocity,(boost::python::arg("idSph")),"Return the normal velocity of the interaction.")
- .def("normalVect",&TemplateFlowEngine::normalVect,(boost::python::arg("idSph")),"Return the normal vector between particles.")
- .def("surfaceDistanceParticle",&TemplateFlowEngine::surfaceDistanceParticle,(boost::python::arg("interaction")),"Return the distance between particles.")
- .def("onlySpheresInteractions",&TemplateFlowEngine::onlySpheresInteractions,(boost::python::arg("interaction")),"Return the id of the interaction only between spheres.")
- .def("pressureProfile",&TemplateFlowEngine::pressureProfile,(boost::python::arg("wallUpY"),boost::python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
- .def("getPorePressure",&TemplateFlowEngine::getPorePressure,(boost::python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
- .def("averageSlicePressure",&TemplateFlowEngine::averageSlicePressure,(boost::python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
- .def("averagePressure",&TemplateFlowEngine::averagePressure,"Measure averaged pore pressure in the entire volume")
- .def("updateBCs",&TemplateFlowEngine::updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
- .def("emulateAction",&TemplateFlowEngine::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
- .def("getCell",&TemplateFlowEngine::getCell,(boost::python::arg("pos")),"get id of the cell containing (X,Y,Z).")
- .def("nCells",&TemplateFlowEngine::nCells,"get the total number of finite cells in the triangulation.")
- .def("getVertices",&TemplateFlowEngine::getVertices,(boost::python::arg("id")),"get the vertices of a cell")
- #ifdef LINSOLV
- .def("exportMatrix",&TemplateFlowEngine::exportMatrix,(boost::python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
- .def("exportTriplets",&TemplateFlowEngine::exportTriplets,(boost::python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
- .def("cholmodStats",&TemplateFlowEngine::cholmodStats,"get statistics of cholmod solver activity")
- .def("metisUsed",&TemplateFlowEngine::metisUsed,"check wether metis lib is effectively used")
- .add_property("forceMetis",&TemplateFlowEngine::getForceMetis,&TemplateFlowEngine::setForceMetis,"If true, METIS is used for matrix preconditioning, else Cholmod is free to choose the best method (which may be METIS to, depending on the matrix). See ``nmethods`` in Cholmod documentation")
- #endif
- .def("compTessVolumes",&TemplateFlowEngine::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
- .def("volume",&TemplateFlowEngine::getVolume,(boost::python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
- .def("averageVelocity",&TemplateFlowEngine::averageVelocity,"measure the mean velocity in the period")
- )
-};
-// Definition of functions in a separate file for clarity
-#include<yade/pkg/pfv/FlowEngine.ipp>
-
-class FlowCellInfo : public CGT::SimpleCellInfo {
- public:
- //For vector storage of all cells
- unsigned int index;
- int volumeSign;
- bool Pcondition;
- Real invVoidV;
- Real t;
- int fict;
- Real volumeVariation;
- double pression;
- //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
- CVector averageCellVelocity;
- // Surface vectors of facets, pointing from outside toward inside the cell
- std::vector<CVector> facetSurfaces;
- //Ratio between fluid surface and facet surface
- std::vector<Real> facetFluidSurfacesRatio;
- // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
- std::vector<CVector> unitForceVectors;
- // Store the area of triangle-sphere intersections for each facet (used in forces definition)
- std::vector<CVector> facetSphereCrossSections;
- std::vector<CVector> cellForce;
- std::vector<double> rayHydr;
- std::vector<double> modulePermeability;
- // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.doInterpolate
- double solidSurfaces [4][4];
-
- FlowCellInfo (void)
- {
- modulePermeability.resize(4, 0);
- cellForce.resize(4,CGAL::NULL_VECTOR);
- facetSurfaces.resize(4,CGAL::NULL_VECTOR);
- facetFluidSurfacesRatio.resize(4,0);
- facetSphereCrossSections.resize(4,CGAL::NULL_VECTOR);
- unitForceVectors.resize(4,CGAL::NULL_VECTOR);
- for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
- rayHydr.resize(4, 0);
- invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;
- isFictious=false; Pcondition = false; isGhost = false;
- isvisited = false;
- isGhost=false;
- }
- bool isGhost;
- double invSumK;
- bool isvisited;
-
- inline Real& volume (void) {return t;}
- inline const Real& invVoidVolume (void) const {return invVoidV;}
- inline Real& invVoidVolume (void) {return invVoidV;}
- inline Real& dv (void) {return volumeVariation;}
- inline int& fictious (void) {return fict;}
- inline double& p (void) {return pression;}
- inline const double shiftedP (void) const {return pression;} //For compatibility with the periodic case
- inline const std::vector<double>& kNorm (void) const {return modulePermeability;}
- inline std::vector<double>& kNorm (void) {return modulePermeability;}
- inline std::vector< CVector >& facetSurf (void) {return facetSurfaces;}
- inline std::vector<CVector>& force (void) {return cellForce;}
- inline std::vector<double>& Rh (void) {return rayHydr;}
- inline CVector& averageVelocity (void) {return averageCellVelocity;}
- //used for transfering values between two triangulations, overload with more variables in derived classes (see e.g. SoluteFlow)
- inline void getInfo(const FlowCellInfo& otherCellInfo) {p()=otherCellInfo.shiftedP();}
-};
-
-class FlowVertexInfo : public CGT::SimpleVertexInfo {
- CVector grainVelocity;
- Real volumeIncidentCells;
-public:
- CVector forces;
- bool isGhost;
- FlowVertexInfo (void) {isGhost=false;}
- inline CVector& force (void) {return forces;}
- inline CVector& vel (void) {return grainVelocity;}
- inline Real& volCells (void) {return volumeIncidentCells;}
- inline const CVector ghostShift (void) const {return CGAL::NULL_VECTOR;}
-};
-
-
-
=== added file 'pkg/pfv/FlowEngine.hpp.in'
--- pkg/pfv/FlowEngine.hpp.in 1970-01-01 00:00:00 +0000
+++ pkg/pfv/FlowEngine.hpp.in 2014-06-20 15:28:21 +0000
@@ -0,0 +1,416 @@
+/*************************************************************************
+* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
+* Copyright (C) 2009 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. *
+*************************************************************************/
+
+/*!
+
+FlowEngine is an interface between Yade and the fluid solvers defined in lib/triangulation and using the PFV scheme.
+There are also a few non-trivial functions defined here, such has building triangulation and computating elements volumes.
+The code strongly relies on CGAL library for building triangulations on the top of sphere packings.
+CGAL's trangulations introduce the data associated to each element of the triangulation through template parameters (Cell_info and Vertex_info).
+FlowEngine is thus defined via the template class TemplateFlowEngine(Cellinfo,VertexInfo), for easier addition of variables in various couplings.
+PeriodicFlowEngine is a variant for periodic boundary conditions.
+
+Which solver will be actually used internally to obtain pore pressure will depend partly on compile time flags (libcholmod available and #define LINSOLV),
+and on runtime settings (useSolver=0: Gauss-Seidel (iterative), useSolver=1: CHOLMOD if available (direct sparse cholesky)).
+
+The files defining lower level classes are in yade/lib. The code uses CGAL::Triangulation3 for managing the mesh and storing data. Eigen3::Sparse, suitesparse::cholmod, and metis are used for solving the linear systems with a direct method (Cholesky). An iterative method (Gauss-Seidel) is implemented directly in Yade and can be used as a fallback (see FlowEngine::useSolver).
+
+Most classes in lib/triangulation are templates, and are therefore completely defined in header files.
+A pseudo hpp/cpp split is reproduced for clarity, with hpp/ipp extensions (but again, in the end they are all include files).
+The files are
+- RegularTriangulation.h : declaration of info types embeded in the mesh (template parameters of CGAL::Triangulation3), and instanciation of RTriangulation classes
+- Tesselation.h/ipp : encapsulate RegularTriangulation and adds functions to manipulate the dual (Voronoi) graph of the triangulation
+- Network.hpp/ipp : specialized for PFV model (the two former are used independently by TesselationWrapper), a set of functions to determine volumes and surfaces of intersections between spheres and various subdomains. Contains two triangulations for smooth transitions while remeshing - e.g. interpolating values in the new mesh using the previous one.
+- FlowBoundingSphere.hpp/ipp and PeriodicFlow.hpp/ipp + LinSolv variants: implement the solver in itself (mesh, boundary conditions, solving, defining fluid-particles interactions)
+- FlowEngine.hpp/ipp/cpp (this file)
+
+Variants for periodic boundary conditions are also present.
+
+*/
+
+#pragma once
+#include<yade/core/PartialEngine.hpp>
+#include<yade/pkg/dem/TriaxialCompressionEngine.hpp>
+#include<yade/pkg/dem/TesselationWrapper.hpp>
+#include<yade/lib/triangulation/FlowBoundingSphere.hpp>
+#include<yade/lib/triangulation/PeriodicFlow.hpp>
+#include<yade/pkg/pfv/FlowEngine.hpp>
+
+template<class _CellInfo, class _VertexInfo, class _Tesselation=CGT::_Tesselation<CGT::TriangulationTypes<_VertexInfo,_CellInfo> >, class solverT=DEFAULTSOLVER >
+class TemplateFlowEngine_@TEMPLATE_FLOW_NAME@ : public PartialEngine
+{
+ public :
+ typedef solverT FlowSolver;
+ typedef FlowSolver Solver;//FIXME: useless alias, search/replace then remove
+ typedef typename FlowSolver::Tesselation Tesselation;
+ typedef typename FlowSolver::RTriangulation RTriangulation;
+ typedef typename FlowSolver::FiniteVerticesIterator FiniteVerticesIterator;
+ typedef typename FlowSolver::FiniteCellsIterator FiniteCellsIterator;
+ typedef typename FlowSolver::CellHandle CellHandle;
+ typedef typename FlowSolver::FiniteEdgesIterator FiniteEdgesIterator;
+ typedef typename FlowSolver::VertexHandle VertexHandle;
+ typedef typename RTriangulation::Triangulation_data_structure::Cell::Info CellInfo;
+ typedef typename RTriangulation::Triangulation_data_structure::Vertex::Info VertexInfo;
+
+// protected:
+ shared_ptr<FlowSolver> solver;
+ shared_ptr<FlowSolver> backgroundSolver;
+ volatile bool backgroundCompleted;
+ Cell cachedCell;
+ struct posData {Body::id_t id; Vector3r pos; Real radius; bool isSphere; bool exists; posData(){exists=0;}};
+ 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
+ virtual void setPositionsBuffer(bool current);
+ virtual void trickPermeability() {};
+
+ public :
+ int retriangulationLastIter;
+ enum {wall_xmin, wall_xmax, wall_ymin, wall_ymax, wall_zmin, wall_zmax};
+ Vector3r normal [6];
+ bool currentTes;
+ bool metisForced;
+ int idOffset;
+ double epsVolCumulative;
+ int ReTrg;
+ int ellapsedIter;
+ void initSolver (FlowSolver& flow);
+ #ifdef LINSOLV
+ void setForceMetis (bool force);
+ bool getForceMetis ();
+ #endif
+ void triangulate (Solver& flow);
+ void addBoundary (Solver& flow);
+ void buildTriangulation (double pZero, Solver& flow);
+ void buildTriangulation (Solver& flow);
+ void updateVolumes (Solver& flow);
+ void initializeVolumes (Solver& flow);
+ void boundaryConditions(Solver& flow);
+ void updateBCs () {
+ if (solver->tesselation().maxId>0) boundaryConditions(*solver);//avoids crash at iteration 0, when the packing is not bounded yet
+ else LOG_ERROR("updateBCs not applied");
+ solver->pressureChanged=true;}
+
+ void imposeFlux(Vector3r pos, Real flux);
+ unsigned int imposePressure(Vector3r pos, Real p);
+ void setImposedPressure(unsigned int cond, Real p);
+ void clearImposedPressure();
+ void clearImposedFlux();
+ void computeViscousForces(Solver& flow);
+ Real getCellFlux(unsigned int cond);
+ Real getBoundaryFlux(unsigned int boundary) {return solver->boundaryFlux(boundary);}
+ Vector3r fluidForce(unsigned int idSph) {
+ const CGT::CVector& f=solver->T[solver->currentTes].vertex(idSph)->info().forces; return Vector3r(f[0],f[1],f[2]);}
+ Vector3r averageVelocity();
+
+ Vector3r shearLubForce(unsigned int id_sph) {
+ return (solver->shearLubricationForces.size()>id_sph)?solver->shearLubricationForces[id_sph]:Vector3r::Zero();}
+ Vector3r shearLubTorque(unsigned int id_sph) {
+ return (solver->shearLubricationTorques.size()>id_sph)?solver->shearLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r pumpLubTorque(unsigned int id_sph) {
+ return (solver->pumpLubricationTorques.size()>id_sph)?solver->pumpLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r twistLubTorque(unsigned int id_sph) {
+ return (solver->twistLubricationTorques.size()>id_sph)?solver->twistLubricationTorques[id_sph]:Vector3r::Zero();}
+ Vector3r normalLubForce(unsigned int id_sph) {
+ return (solver->normalLubricationForce.size()>id_sph)?solver->normalLubricationForce[id_sph]:Vector3r::Zero();}
+ Matrix3r bodyShearLubStress(unsigned int id_sph) {
+ return (solver->shearLubricationBodyStress.size()>id_sph)?solver->shearLubricationBodyStress[id_sph]:Matrix3r::Zero();}
+ Matrix3r bodyNormalLubStress(unsigned int id_sph) {
+ return (solver->normalLubricationBodyStress.size()>id_sph)?solver->normalLubricationBodyStress[id_sph]:Matrix3r::Zero();}
+ Vector3r shearVelocity(unsigned int interaction) {
+ return (solver->deltaShearVel[interaction]);}
+ Vector3r normalVelocity(unsigned int interaction) {
+ return (solver->deltaNormVel[interaction]);}
+ Matrix3r normalStressInteraction(unsigned int interaction) {
+ return (solver->normalStressInteraction[interaction]);}
+ Matrix3r shearStressInteraction(unsigned int interaction) {
+ return (solver->shearStressInteraction[interaction]);}
+ Vector3r normalVect(unsigned int interaction) {
+ return (solver->normalV[interaction]);}
+ Real surfaceDistanceParticle(unsigned int interaction) {
+ return (solver->surfaceDistance[interaction]);}
+ Real edgeSize() {
+ return (solver->edgeIds.size());}
+ Real OSI() {
+ return (solver->onlySpheresInteractions.size());}
+ int onlySpheresInteractions(unsigned int interaction) {
+ return (solver->onlySpheresInteractions[interaction]);}
+ boost::python::list getConstrictions(bool all) {
+ vector<Real> csd=solver->getConstrictions(); boost::python::list pycsd;
+ for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k]!=0) || csd[k]>0) pycsd.append(csd[k]); return pycsd;}
+ boost::python::list getConstrictionsFull(bool all) {
+ vector<Constriction> csd=solver->getConstrictionsFull(); boost::python::list pycsd;
+ for (unsigned int k=0;k<csd.size();k++) if ((all && csd[k].second[0]!=0) || csd[k].second[0]>0) {
+ boost::python::list cons;
+ cons.append(csd[k].first.first);
+ cons.append(csd[k].first.second);
+ cons.append(csd[k].second[0]);
+ cons.append(csd[k].second[1]);
+ cons.append(csd[k].second[2]);
+ cons.append(csd[k].second[3]);
+ pycsd.append(cons);}
+ return pycsd;}
+
+ template<class Cellhandle>
+ Real volumeCellSingleFictious (Cellhandle cell);
+ template<class Cellhandle>
+ Real volumeCellDoubleFictious (Cellhandle cell);
+ template<class Cellhandle>
+ Real volumeCellTripleFictious (Cellhandle cell);
+ template<class Cellhandle>
+ Real volumeCell (Cellhandle cell);
+ void Oedometer_Boundary_Conditions();
+ void averageRealCellVelocity();
+ void saveVtk(const char* folder) {solver->saveVtk(folder);}
+ vector<Real> avFlVelOnSph(unsigned int idSph) {return solver->averageFluidVelocityOnSphere(idSph);}
+
+// void setBoundaryVel(Vector3r vel) {topBoundaryVelocity=vel; updateTriangulation=true;}
+ void pressureProfile(double wallUpY, double wallDownY) {return solver->measurePressureProfile(wallUpY,wallDownY);}
+ double getPorePressure(Vector3r pos){return solver->getPorePressure(pos[0], pos[1], pos[2]);}
+ int getCell(double posX, double posY, double posZ){return solver->getCell(posX, posY, posZ);}
+ unsigned int nCells(){return solver->T[solver->currentTes].cellHandles.size();}
+ boost::python::list getVertices(unsigned int id){
+ boost::python::list ids;
+ if (id>=solver->T[solver->currentTes].cellHandles.size()) {LOG_ERROR("id out of range, max value is "<<solver->T[solver->currentTes].cellHandles.size()); return ids;}
+ for (unsigned int i=0;i<4;i++) ids.append(solver->T[solver->currentTes].cellHandles[id]->vertex(i)->info().id());
+ return ids;
+ }
+ double averageSlicePressure(double posY){return solver->averageSlicePressure(posY);}
+ double averagePressure(){return solver->averagePressure();}
+
+ void emulateAction(){
+ scene = Omega::instance().getScene().get();
+ action();}
+ #ifdef LINSOLV
+ void exportMatrix(string filename) {if (useSolver==3) solver->exportMatrix(filename.c_str());
+ else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
+ void exportTriplets(string filename) {if (useSolver==3) solver->exportTriplets(filename.c_str());
+ else cerr<<"available for Cholmod solver (useSolver==3)"<<endl;}
+ void cholmodStats() {
+ cerr << cholmod_print_common((char*)string("PFV Cholmod factorization").c_str(),&(solver->eSolver.cholmod()))<<endl;
+ cerr << "cholmod method:" << solver->eSolver.cholmod().selected<<endl;
+ cerr << "METIS called:"<<solver->eSolver.cholmod().called_nd<<endl;}
+ bool metisUsed() {return bool(solver->eSolver.cholmod().called_nd);}
+ #endif
+
+ virtual ~TemplateFlowEngine_@TEMPLATE_FLOW_NAME@();
+ virtual void action();
+ virtual void backgroundAction();
+
+ //commodities
+ void compTessVolumes() {
+ solver->T[solver->currentTes].compute();
+ solver->T[solver->currentTes].computeVolumes();
+ }
+ Real getVolume (Body::id_t id) {
+ if (solver->T[solver->currentTes].Max_id() <= 0) {emulateAction(); /*LOG_WARN("Not triangulated yet, emulating action");*/}
+ if (solver->T[solver->currentTes].Volume(id) == -1) {compTessVolumes();/* LOG_WARN("Computing all volumes now, as you did not request it explicitely.");*/}
+ return (solver->T[solver->currentTes].Max_id() >= id) ? solver->T[solver->currentTes].Volume(id) : -1;}
+
+ YADE_CLASS_PYCLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(TemplateFlowEngine_@TEMPLATE_FLOW_NAME@,@TEMPLATE_FLOW_NAME@,PartialEngine,"An engine to solve flow problem in saturated granular media. Model description can be found in [Chareyre2012a]_ and [Catalano2014a]_. See the example script FluidCouplingPFV/oedometer.py. More documentation to come.\n\n.. note::Multi-threading seems to work fine for Cholesky decomposition, but it fails for the solve phase in which -j1 is the fastest, here we specify thread numbers independently using :yref:`FlowEngine::numFactorizeThreads` and :yref:`FlowEngine::numSolveThreads`. These multhreading settings are only impacting the behaviour of openblas library and are relatively independant of :yref:`FlowEngine::multithread`. However, the settings have to be globally consistent. For instance, :yref:`multithread<FlowEngine::multithread>` =True with yref:`numFactorizeThreads<FlowEngine::numFactorizeThreads>` = yref:`numSolveThreads<FlowEngine::numSolveThreads>` = 4 implies that openblas will mobilize 8 processors at some point. If the system does not have so many procs. it will hurt performance.",
+ ((bool,isActivated,true,,"Activates Flow Engine"))
+ ((bool,first,true,,"Controls the initialization/update phases"))
+ ((bool,doInterpolate,false,,"Force the interpolation of cell's info while remeshing. By default, interpolation would be done only for compressible fluids. It can be forced with this flag."))
+ ((double, fluidBulkModulus, 0.,,"Bulk modulus of fluid (inverse of compressibility) K=-dP*V/dV [Pa]. Flow is compressible if fluidBulkModulus > 0, else incompressible."))
+ ((Real, dt, 0,,"timestep [s]"))
+ ((bool,permeabilityMap,false,,"Enable/disable stocking of average permeability scalar in cell infos."))
+ ((bool, slipBoundary, true,, "Controls friction condition on lateral walls"))
+ ((bool,waveAction, false,, "Allow sinusoidal pressure condition to simulate ocean waves"))
+ ((double, sineMagnitude, 0,, "Pressure value (amplitude) when sinusoidal pressure is applied (p )"))
+ ((double, sineAverage, 0,,"Pressure value (average) when sinusoidal pressure is applied"))
+ ((bool, debug, false,,"Activate debug messages"))
+ ((double, wallThickness,0.001,,"Walls thickness"))
+ ((double,pZero,0,,"The value used for initializing pore pressure. It is useless for incompressible fluid, but important for compressible model."))
+ ((double,tolerance,1e-06,,"Gauss-Seidel tolerance"))
+ ((double,relax,1.9,,"Gauss-Seidel relaxation"))
+ ((bool, updateTriangulation, 0,,"If true the medium is retriangulated. Can be switched on to force retriangulation after some events (else it will be true periodicaly based on :yref:`FlowEngine::defTolerance` and :yref:`FlowEngine::meshUpdateInterval`. Of course, it costs CPU time."))
+ ((int,meshUpdateInterval,1000,,"Maximum number of timesteps between re-triangulation events. See also :yref:`FlowEngine::defTolerance`."))
+ ((double, epsVolMax, 0,(Attr::readonly),"Maximal absolute volumetric strain computed at each iteration. |yupdate|"))
+ ((double, defTolerance,0.05,,"Cumulated deformation threshold for which retriangulation of pore space is performed. If negative, the triangulation update will occure with a fixed frequency on the basis of :yref:`FlowEngine::meshUpdateInterval`"))
+ ((double, porosity, 0,(Attr::readonly),"Porosity computed at each retriangulation |yupdate|"))
+ ((bool,meanKStat,false,,"report the local permeabilities' correction"))
+ ((bool,clampKValues,true,,"If true, clamp local permeabilities in [minKdivKmean,maxKdivKmean]*globalK. This clamping can avoid singular values in the permeability matrix and may reduce numerical errors in the solve phase. It will also hide junk values if they exist, or bias all values in very heterogeneous problems. So, use this with care."))
+ ((Real,minKdivKmean,0.0001,,"define the min K value (see :yref:`FlowEngine::clampKValues`)"))
+ ((Real,maxKdivKmean,100,,"define the max K value (see :yref:`FlowEngine::clampKValues`)"))
+ ((double,permeabilityFactor,1.0,,"permability multiplier"))
+ ((double,viscosity,1.0,,"viscosity of the fluid"))
+ ((double,stiffness, 10000,,"equivalent contact stiffness used in the lubrication model"))
+ ((int, useSolver, 0,, "Solver to use 0=G-Seidel, 1=Taucs, 2-Pardiso, 3-CHOLMOD"))
+ ((int, xmin,0,(Attr::readonly),"Index of the boundary $x_{min}$. This index is not equal the the id of the corresponding body in general, it may be used to access the corresponding attributes (e.g. flow.bndCondValue[flow.xmin], flow.wallId[flow.xmin],...)."))
+ ((int, xmax,1,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, ymin,2,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, ymax,3,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, zmin,4,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+ ((int, zmax,5,(Attr::readonly),"See :yref:`FlowEngine::xmin`."))
+
+ ((vector<bool>, bndCondIsPressure, vector<bool>(6,false),,"defines the type of boundary condition for each side. True if pressure is imposed, False for no-flux. Indexes can be retrieved with :yref:`FlowEngine::xmin` and friends."))
+ ((vector<double>, bndCondValue, vector<double>(6,0),,"Imposed value of a boundary condition. Only applies if the boundary condition is imposed pressure, else the imposed flux is always zero presently (may be generalized to non-zero imposed fluxes in the future)."))
+ //FIXME: to be implemented:
+ ((vector<Vector3r>, boundaryVelocity, vector<Vector3r>(6,Vector3r::Zero()),, "velocity on top boundary, only change it using :yref:`FlowEngine::setBoundaryVel`"))
+ ((int, ignoredBody,-1,,"Id of a sphere to exclude from the triangulation.)"))
+ ((vector<int>, wallIds,vector<int>(6),,"body ids of the boundaries (default values are ok only if aabbWalls are appended before spheres, i.e. numbered 0,...,5)"))
+ ((vector<bool>, boundaryUseMaxMin, vector<bool>(6,true),,"If true (default value) bounding sphere is added as function of max/min sphere coord, if false as function of yade wall position"))
+// ((bool, display_force, false,,"display the lubrication force applied on particles"))
+// ((bool, create_file, false,,"create file of velocities"))
+ ((bool, viscousShear, false,,"compute viscous shear terms as developped by Donia Marzougui (FIXME: ref.)"))
+ ((bool, shearLubrication, false,,"compute shear lubrication force as developped by Brule (FIXME: ref.) "))
+ ((bool, pumpTorque, false,,"Compute pump torque applied on particles "))
+ ((bool, twistTorque, false,,"Compute twist torque applied on particles "))
+ ((double, eps, 0.00001,,"roughness defined as a fraction of particles size, giving the minimum distance between particles in the lubrication model."))
+ ((bool, pressureForce, true,,"compute the pressure field and associated fluid forces. WARNING: turning off means fluid flow is not computed at all."))
+ ((bool, normalLubrication, false,,"compute normal lubrication force as developped by Brule"))
+ ((bool, viscousNormalBodyStress, false,,"compute normal viscous stress applied on each body"))
+ ((bool, viscousShearBodyStress, false,,"compute shear viscous stress applied on each body"))
+ ((bool, multithread, false,,"Build triangulation and factorize in the background (multi-thread mode)"))
+ #ifdef EIGENSPARSE_LIB
+ ((int, numSolveThreads, 1,,"number of openblas threads in the solve phase."))
+ ((int, numFactorizeThreads, 1,,"number of openblas threads in the factorization phase"))
+ #endif
+ ((vector<Real>, boundaryPressure,vector<Real>(),,"values defining pressure along x-axis for the top surface. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryXPos`"))
+ ((vector<Real>, boundaryXPos,vector<Real>(),,"values of the x-coordinate for which pressure is defined. See also :yref:`@TEMPLATE_FLOW_NAME@::boundaryPressure`"))
+ ,
+ /*deprec*/
+ ((meanK_opt,clampKValues,"the name changed"))
+ ,,
+ timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
+ for (int i=0; i<6; ++i){normal[i]=Vector3r::Zero(); wallIds[i]=i;}
+ normal[wall_ymin].y()=normal[wall_xmin].x()=normal[wall_zmin].z()=1;
+ normal[wall_ymax].y()=normal[wall_xmax].x()=normal[wall_zmax].z()=-1;
+ solver = shared_ptr<FlowSolver> (new FlowSolver);
+ first=true;
+ epsVolMax=epsVolCumulative=retriangulationLastIter=0;
+ ReTrg=1;
+ backgroundCompleted=true;
+ ellapsedIter=0;
+ metisForced=false;
+ ,
+ .def("imposeFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::imposeFlux,(boost::python::arg("pos"),boost::python::arg("p")),"Impose incoming flux in boundary cell of location 'pos'.")
+ .def("imposePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::imposePressure,(boost::python::arg("pos"),boost::python::arg("p")),"Impose pressure in cell of location 'pos'. The index of the condition is returned (for multiple imposed pressures at different points).")
+ .def("setImposedPressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::setImposedPressure,(boost::python::arg("cond"),boost::python::arg("p")),"Set pressure value at the point indexed 'cond'.")
+ .def("clearImposedPressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::clearImposedPressure,"Clear the list of points with pressure imposed.")
+ .def("clearImposedFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::clearImposedFlux,"Clear the list of points with flux imposed.")
+ .def("getCellFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCellFlux,(boost::python::arg("cond")),"Get influx in cell associated to an imposed P (indexed using 'cond').")
+ .def("getBoundaryFlux",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getBoundaryFlux,(boost::python::arg("boundary")),"Get total flux through boundary defined by its body id.\n\n.. note:: The flux may be not zero even for no-flow condition. This artifact comes from cells which are incident to two or more boundaries (along the edges of the sample, typically). Such flux evaluation on impermeable boundary is just irrelevant, it does not imply that the boundary condition is not applied properly.")
+ .def("getConstrictions",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getConstrictions,(boost::python::arg("all")=true),"Get the list of constriction radii (inscribed circle) for all finite facets (if all==True) or all facets not incident to a virtual bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres.")
+ .def("getConstrictionsFull",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getConstrictionsFull,(boost::python::arg("all")=true),"Get the list of constrictions (inscribed circle) for all finite facets (if all==True), or all facets not incident to a fictious bounding sphere (if all==False). When all facets are returned, negative radii denote facet incident to one or more fictious spheres. The constrictions are returned in the format {{cell1,cell2}{rad,nx,ny,nz}}")
+ .def("edgeSize",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::edgeSize,"Return the number of interactions.")
+ .def("OSI",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::OSI,"Return the number of interactions only between spheres.")
+
+ .def("saveVtk",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::saveVtk,(boost::python::arg("folder")="./VTK"),"Save pressure field in vtk format. Specify a folder name for output.")
+ .def("avFlVelOnSph",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::avFlVelOnSph,(boost::python::arg("idSph")),"compute a sphere-centered average fluid velocity")
+ .def("fluidForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::fluidForce,(boost::python::arg("idSph")),"Return the fluid force on sphere idSph.")
+ .def("shearLubForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::shearLubForce,(boost::python::arg("idSph")),"Return the shear lubrication force on sphere idSph.")
+ .def("shearLubTorque",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::shearLubTorque,(boost::python::arg("idSph")),"Return the shear lubrication torque on sphere idSph.")
+ .def("normalLubForce",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::normalLubForce,(boost::python::arg("idSph")),"Return the normal lubrication force on sphere idSph.")
+ .def("bodyShearLubStress",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::bodyShearLubStress,(boost::python::arg("idSph")),"Return the shear lubrication stress on sphere idSph.")
+ .def("bodyNormalLubStress",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::bodyNormalLubStress,(boost::python::arg("idSph")),"Return the normal lubrication stress on sphere idSph.")
+ .def("shearVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::shearVelocity,(boost::python::arg("idSph")),"Return the shear velocity of the interaction.")
+ .def("normalVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::normalVelocity,(boost::python::arg("idSph")),"Return the normal velocity of the interaction.")
+ .def("normalVect",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::normalVect,(boost::python::arg("idSph")),"Return the normal vector between particles.")
+ .def("surfaceDistanceParticle",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::surfaceDistanceParticle,(boost::python::arg("interaction")),"Return the distance between particles.")
+ .def("onlySpheresInteractions",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::onlySpheresInteractions,(boost::python::arg("interaction")),"Return the id of the interaction only between spheres.")
+ .def("pressureProfile",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::pressureProfile,(boost::python::arg("wallUpY"),boost::python::arg("wallDownY")),"Measure pore pressure in 6 equally-spaced points along the height of the sample")
+ .def("getPorePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getPorePressure,(boost::python::arg("pos")),"Measure pore pressure in position pos[0],pos[1],pos[2]")
+ .def("averageSlicePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageSlicePressure,(boost::python::arg("posY")),"Measure slice-averaged pore pressure at height posY")
+ .def("averagePressure",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averagePressure,"Measure averaged pore pressure in the entire volume")
+ .def("updateBCs",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::updateBCs,"tells the engine to update it's boundary conditions before running (especially useful when changing boundary pressure - should not be needed for point-wise imposed pressure)")
+ .def("emulateAction",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::emulateAction,"get scene and run action (may be used to manipulate an engine outside the timestepping loop).")
+ .def("getCell",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getCell,(boost::python::arg("pos")),"get id of the cell containing (X,Y,Z).")
+ .def("nCells",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::nCells,"get the total number of finite cells in the triangulation.")
+ .def("getVertices",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVertices,(boost::python::arg("id")),"get the vertices of a cell")
+ #ifdef LINSOLV
+ .def("exportMatrix",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::exportMatrix,(boost::python::arg("filename")="matrix"),"Export system matrix to a file with all entries (even zeros will displayed).")
+ .def("exportTriplets",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::exportTriplets,(boost::python::arg("filename")="triplets"),"Export system matrix to a file with only non-zero entries.")
+ .def("cholmodStats",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::cholmodStats,"get statistics of cholmod solver activity")
+ .def("metisUsed",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::metisUsed,"check wether metis lib is effectively used")
+ .add_property("forceMetis",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getForceMetis,&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::setForceMetis,"If true, METIS is used for matrix preconditioning, else Cholmod is free to choose the best method (which may be METIS to, depending on the matrix). See ``nmethods`` in Cholmod documentation")
+ #endif
+ .def("compTessVolumes",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::compTessVolumes,"Like TesselationWrapper::computeVolumes()")
+ .def("volume",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::getVolume,(boost::python::arg("id")=0),"Returns the volume of Voronoi's cell of a sphere.")
+ .def("averageVelocity",&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::averageVelocity,"measure the mean velocity in the period")
+ )
+};
+// Definition of functions in a separate file for clarity
+
+class FlowCellInfo_@TEMPLATE_FLOW_NAME@ : public CGT::SimpleCellInfo {
+ public:
+ //For vector storage of all cells
+ unsigned int index;
+ int volumeSign;
+ bool Pcondition;
+ Real invVoidV;
+ Real t;
+ int fict;
+ Real volumeVariation;
+ double pression;
+ //average relative (fluid - facet translation) velocity defined for a single cell as 1/Volume * SUM_ON_FACETS(x_average_facet*average_facet_flow_rate)
+ CVector averageCellVelocity;
+ // Surface vectors of facets, pointing from outside toward inside the cell
+ std::vector<CVector> facetSurfaces;
+ //Ratio between fluid surface and facet surface
+ std::vector<Real> facetFluidSurfacesRatio;
+ // Reflects the geometrical property of the cell, so that the force by cell fluid on grain "i" is pressure*unitForceVectors[i]
+ std::vector<CVector> unitForceVectors;
+ // Store the area of triangle-sphere intersections for each facet (used in forces definition)
+ std::vector<CVector> facetSphereCrossSections;
+ std::vector<CVector> cellForce;
+ std::vector<double> rayHydr;
+ std::vector<double> modulePermeability;
+ // Partial surfaces of spheres in the double-tetrahedron linking two voronoi centers. [i][j] is for sphere facet "i" and sphere facetVertices[i][j]. Last component for 1/sum_surfaces in the facet.doInterpolate
+ double solidSurfaces [4][4];
+
+ FlowCellInfo_@TEMPLATE_FLOW_NAME@ (void)
+ {
+ modulePermeability.resize(4, 0);
+ cellForce.resize(4,CGAL::NULL_VECTOR);
+ facetSurfaces.resize(4,CGAL::NULL_VECTOR);
+ facetFluidSurfacesRatio.resize(4,0);
+ facetSphereCrossSections.resize(4,CGAL::NULL_VECTOR);
+ unitForceVectors.resize(4,CGAL::NULL_VECTOR);
+ for (int k=0; k<4;k++) for (int l=0; l<3;l++) solidSurfaces[k][l]=0;
+ rayHydr.resize(4, 0);
+ invSumK=index=volumeSign=s=volumeVariation=pression=invVoidV=fict=0;
+ isFictious=false; Pcondition = false; isGhost = false;
+ isvisited = false;
+ isGhost=false;
+ }
+ bool isGhost;
+ double invSumK;
+ bool isvisited;
+
+ inline Real& volume (void) {return t;}
+ inline const Real& invVoidVolume (void) const {return invVoidV;}
+ inline Real& invVoidVolume (void) {return invVoidV;}
+ inline Real& dv (void) {return volumeVariation;}
+ inline int& fictious (void) {return fict;}
+ inline double& p (void) {return pression;}
+ inline const double shiftedP (void) const {return pression;} //For compatibility with the periodic case
+ inline const std::vector<double>& kNorm (void) const {return modulePermeability;}
+ inline std::vector<double>& kNorm (void) {return modulePermeability;}
+ inline std::vector< CVector >& facetSurf (void) {return facetSurfaces;}
+ inline std::vector<CVector>& force (void) {return cellForce;}
+ inline std::vector<double>& Rh (void) {return rayHydr;}
+ inline CVector& averageVelocity (void) {return averageCellVelocity;}
+ //used for transfering values between two triangulations, overload with more variables in derived classes (see e.g. SoluteFlow)
+ inline void getInfo(const FlowCellInfo_@TEMPLATE_FLOW_NAME@& otherCellInfo) {p()=otherCellInfo.shiftedP();}
+};
+
+class FlowVertexInfo_@TEMPLATE_FLOW_NAME@ : public CGT::SimpleVertexInfo {
+ CVector grainVelocity;
+ Real volumeIncidentCells;
+public:
+ CVector forces;
+ bool isGhost;
+ FlowVertexInfo_@TEMPLATE_FLOW_NAME@ (void) {isGhost=false;}
+ inline CVector& force (void) {return forces;}
+ inline CVector& vel (void) {return grainVelocity;}
+ inline Real& volCells (void) {return volumeIncidentCells;}
+ inline const CVector ghostShift (void) const {return CGAL::NULL_VECTOR;}
+};
+
+#include "FlowEngine_@TEMPLATE_FLOW_NAME@.ipp"
=== modified file 'pkg/pfv/FlowEngine.ipp'
--- pkg/pfv/FlowEngine.ipp 2014-06-04 17:19:50 +0000
+++ pkg/pfv/FlowEngine.ipp 2014-06-20 15:28:21 +0000
@@ -1,750 +1,1 @@
-/*************************************************************************
-* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
-* Copyright (C) 2009 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. *
-*************************************************************************/
-#ifdef YADE_CGAL
-
-#ifdef FLOW_ENGINE
-#include<yade/core/Scene.hpp>
-#include<yade/lib/base/Math.hpp>
-#include<yade/pkg/dem/TesselationWrapper.hpp>
-#include<yade/pkg/common/Sphere.hpp>
-#include<yade/pkg/common/Wall.hpp>
-#include<yade/pkg/common/Box.hpp>
-#include <sys/stat.h>
-#include <sys/types.h>
-#include <boost/thread.hpp>
-#include <boost/date_time.hpp>
-#include <boost/date_time/posix_time/posix_time.hpp>
-
-#ifdef LINSOLV
-#include <cholmod.h>
-#endif
-
-#include "FlowEngine.hpp"
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine() {}
-
-// YADE_PLUGIN((TFlowEng));
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-unsigned int TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposePressure(Vector3r pos, Real p)
-{
-// if (!flow) LOG_ERROR("no flow defined yet, run at least one iter");
- solver->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) );
- //force immediate update of boundary conditions
- updateTriangulation=true;
- return solver->imposedP.size()-1;
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::action()
-{
- if ( !isActivated ) return;
- timingDeltas->start();
- setPositionsBuffer(true);
- timingDeltas->checkpoint ( "Position buffer" );
- if (first) {
- if (multithread) setPositionsBuffer(false);
- buildTriangulation(pZero,*solver);
- initializeVolumes(*solver);
- backgroundSolver=solver;
- backgroundCompleted=true;
- }
- solver->ompThreads = ompThreads>0? ompThreads : omp_get_max_threads();
-
- timingDeltas->checkpoint ( "Triangulating" );
- updateVolumes ( *solver );
- timingDeltas->checkpoint ( "Update_Volumes" );
-
- epsVolCumulative += epsVolMax;
- retriangulationLastIter++;
- if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
- (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
-
- ///compute flow and and forces here
- if (pressureForce){
- solver->gaussSeidel(scene->dt);
- timingDeltas->checkpoint ( "Gauss-Seidel (includes matrix construct and factorization in single-thread mode)" );
- solver->computeFacetForcesWithCache();}
- timingDeltas->checkpoint ( "compute_Forces" );
- ///Application of vicscous forces
- scene->forces.sync();
- timingDeltas->checkpoint ( "forces.sync()" );
- computeViscousForces ( *solver );
- timingDeltas->checkpoint ( "viscous forces" );
- Vector3r force;
- Vector3r torque;
- FiniteVerticesIterator verticesEnd = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
- for ( FiniteVerticesIterator vIt = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); vIt != verticesEnd; vIt++ ) {
- force = pressureForce ? Vector3r ( vIt->info().forces[0],vIt->info().forces[1],vIt->info().forces[2] ): Vector3r(0,0,0);
- torque = Vector3r(0,0,0);
- if (shearLubrication || viscousShear){
- force = force + solver->shearLubricationForces[vIt->info().id()];
- torque = torque + solver->shearLubricationTorques[vIt->info().id()];
- if (pumpTorque)
- torque = torque + solver->pumpLubricationTorques[vIt->info().id()];
- }
- if (twistTorque)
- torque = torque + solver->twistLubricationTorques[vIt->info().id()];
- if (normalLubrication)
- force = force + solver-> normalLubricationForce[vIt->info().id()];
- scene->forces.addForce ( vIt->info().id(), force);
- scene->forces.addTorque ( vIt->info().id(), torque);
- }
- ///End compute flow and forces
- timingDeltas->checkpoint ( "Applying Forces" );
- #ifdef LINSOLV
- int sleeping = 0;
- if (multithread && !first) {
- while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/
- sleeping++;
- boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
- if (debug && sleeping) cerr<<"sleeping..."<<sleeping<<endl;
- if (updateTriangulation || (ellapsedIter>(0.5*meshUpdateInterval) && backgroundCompleted)) {
- if (debug) cerr<<"switch flow solver"<<endl;
- if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
- if (fluidBulkModulus>0 || doInterpolate) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
- solver=backgroundSolver;
- backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
- if (metisForced) {backgroundSolver->eSolver.cholmod().nmethods=1; backgroundSolver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;}
- //Copy imposed pressures/flow from the old solver
- backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
- backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
- if (debug) cerr<<"switched"<<endl;
- setPositionsBuffer(false);//set "parallel" buffer for background calculation
- backgroundCompleted=false;
- retriangulationLastIter=ellapsedIter;
- updateTriangulation=false;
- epsVolCumulative=0;
- ellapsedIter=0;
- boost::thread workerThread(&TemplateFlowEngine::backgroundAction,this);
- workerThread.detach();
- if (debug) cerr<<"backgrounded"<<endl;
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- if (debug) cerr<<"volumes initialized"<<endl;
- }
- else {
- if (debug && !backgroundCompleted) cerr<<"still computing solver in the background, ellapsedIter="<<ellapsedIter<<endl;
- ellapsedIter++;
- }
- } else
- #endif
- {
- if (updateTriangulation && !first) {
- buildTriangulation (pZero, *solver);
- initializeVolumes(*solver);
- computeViscousForces(*solver);
- updateTriangulation = false;
- epsVolCumulative=0;
- retriangulationLastIter=0;
- ReTrg++;}
- }
- first=false;
- timingDeltas->checkpoint ( "triangulate + init volumes" );
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::backgroundAction()
-{
- if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
- buildTriangulation ( pZero,*backgroundSolver );
- //FIXME: GS is computing too much, we need only matrix factorization in fact
- backgroundSolver->gaussSeidel(scene->dt);
- //FIXME(2): and here we need only cached variables, not forces
- backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
-// boost::this_thread::sleep(boost::posix_time::seconds(5));
- backgroundCompleted = true;
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::boundaryConditions ( Solver& flow )
-{
- for (int k=0;k<6;k++) {
- flow.boundary (wallIds[k]).flowCondition=!bndCondIsPressure[k];
- flow.boundary (wallIds[k]).value=bndCondValue[k];
- flow.boundary (wallIds[k]).velocity = boundaryVelocity[k];//FIXME: needs correct implementation, maybe update the cached pos/vel?
- }
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setImposedPressure ( unsigned int cond, Real p)
-{
- if ( cond>=solver->imposedP.size() ) LOG_ERROR ( "Setting p with cond higher than imposedP size." );
- solver->imposedP[cond].second=p;
- //force immediate update of boundary conditions
- solver->pressureChanged=true;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposeFlux ( Vector3r pos, Real flux){
- solver->imposedF.push_back ( pair<CGT::Point,Real> ( CGT::Point ( pos[0],pos[1],pos[2] ), flux ) );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedPressure () { solver->imposedP.clear(); solver->IPCells.clear();}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedFlux () { solver->imposedF.clear(); solver->IFCells.clear();}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getCellFlux ( unsigned int cond)
-{
- if ( cond>=solver->imposedP.size() ) {LOG_ERROR ( "Getting flux with cond higher than imposedP size." ); return 0;}
- double flux=0;
- typename Solver::CellHandle& cell= solver->IPCells[cond];
- for ( int ngb=0;ngb<4;ngb++ ) {
- /*if (!cell->neighbor(ngb)->info().Pcondition)*/
- flux+= cell->info().kNorm() [ngb]* ( cell->info().p()-cell->neighbor ( ngb )->info().p() );
- }
- return flux+cell->info().dv();
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initSolver ( FlowSolver& flow )
-{
- flow.Vtotalissimo=0; flow.VSolidTot=0; flow.vPoral=0; flow.sSolidTot=0;
- flow.slipBoundary=slipBoundary;
- flow.kFactor = permeabilityFactor;
- flow.debugOut = debug;
- flow.useSolver = useSolver;
- #ifdef EIGENSPARSE_LIB
- flow.numSolveThreads = numSolveThreads;
- flow.numFactorizeThreads = numFactorizeThreads;
- #endif
- flow.meanKStat = meanKStat;
- flow.viscosity = viscosity;
- flow.tolerance=tolerance;
- flow.relax=relax;
- flow.clampKValues = clampKValues;
- flow.maxKdivKmean = maxKdivKmean;
- flow.minKdivKmean = minKdivKmean;
- flow.meanKStat = meanKStat;
- flow.permeabilityMap = permeabilityMap;
- flow.fluidBulkModulus = fluidBulkModulus;
-// flow.tesselation().Clear();
- flow.tesselation().maxId=-1;
- flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
-}
-
-#ifdef LINSOLV
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setForceMetis ( bool force )
-{
- if (force) {
- metisForced=true;
- solver->eSolver.cholmod().nmethods=1;
- solver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;
- } else {cholmod_defaults(&(solver->eSolver.cholmod())); metisForced=false;}
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-bool TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::getForceMetis () {return (solver->eSolver.cholmod().nmethods==1);}
-#endif
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( Solver& flow )
-{
- buildTriangulation ( 0.f,flow );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-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;}
- flow.resetNetwork();
- initSolver(flow);
-
- addBoundary ( flow );
- triangulate ( flow );
- if ( debug ) cout << endl << "Tesselating------" << endl << endl;
- flow.tesselation().compute();
-
- flow.defineFictiousCells();
- // For faster loops on cells define this vector
- flow.tesselation().cellHandles.clear();
- flow.tesselation().cellHandles.reserve(flow.tesselation().Triangulation().number_of_finite_cells());
- FiniteCellsIterator cell_end = flow.tesselation().Triangulation().finite_cells_end();
- int k=0;
- for ( FiniteCellsIterator cell = flow.tesselation().Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
- flow.tesselation().cellHandles.push_back(cell);
- cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
- flow.displayStatistics ();
- flow.computePermeability();
- //This virtual function does nothing yet, derived class may overload it to make permeability different (see DFN engine)
- trickPermeability();
- porosity = flow.vPoralPorosity/flow.vTotalPorosity;
-
- boundaryConditions ( flow );
- flow.initializePressure ( pZero );
-
- if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.tesselation() );
- if ( waveAction ) flow.applySinusoidalPressure ( flow.tesselation().Triangulation(), sineMagnitude, sineAverage, 30 );
- else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.tesselation().Triangulation(), boundaryXPos , boundaryPressure);
- if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::setPositionsBuffer(bool current)
-{
- vector<posData>& buffer = current? positionBufferCurrent : positionBufferParallel;
- 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;
- }
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::addBoundary ( Solver& flow )
-{
- vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
- solver->xMin = Mathr::MAX_REAL, solver->xMax = -Mathr::MAX_REAL, solver->yMin = Mathr::MAX_REAL, solver->yMax = -Mathr::MAX_REAL, solver->zMin = Mathr::MAX_REAL, solver->zMax = -Mathr::MAX_REAL;
- FOREACH ( const posData& b, buffer ) {
- if ( !b.exists ) continue;
- if ( b.isSphere ) {
- const Real& rad = b.radius;
- const Real& x = b.pos[0];
- const Real& y = b.pos[1];
- const Real& z = b.pos[2];
- flow.xMin = min ( flow.xMin, x-rad );
- flow.xMax = max ( flow.xMax, x+rad );
- flow.yMin = min ( flow.yMin, y-rad );
- flow.yMax = max ( flow.yMax, y+rad );
- flow.zMin = min ( flow.zMin, z-rad );
- flow.zMax = max ( flow.zMax, z+rad );
- }
- }
- //FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
- idOffset = flow.tesselation().maxId+1;
- flow.idOffset = idOffset;
- flow.sectionArea = ( flow.xMax - flow.xMin ) * ( flow.zMax-flow.zMin );
- flow.vTotal = ( flow.xMax-flow.xMin ) * ( flow.yMax-flow.yMin ) * ( flow.zMax-flow.zMin );
- flow.yMinId=wallIds[ymin];
- flow.yMaxId=wallIds[ymax];
- flow.xMaxId=wallIds[xmax];
- flow.xMinId=wallIds[xmin];
- flow.zMinId=wallIds[zmin];
- flow.zMaxId=wallIds[zmax];
-
- //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
- flow.boundsIds[0]= &flow.xMinId;
- flow.boundsIds[1]= &flow.xMaxId;
- flow.boundsIds[2]= &flow.yMinId;
- flow.boundsIds[3]= &flow.yMaxId;
- flow.boundsIds[4]= &flow.zMinId;
- flow.boundsIds[5]= &flow.zMaxId;
-
- for (int k=0;k<6;k++) flow.boundary ( *flow.boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k];
-
- flow.cornerMin = CGT::Point ( flow.xMin, flow.yMin, flow.zMin );
- flow.cornerMax = CGT::Point ( flow.xMax, flow.yMax, flow.zMax );
-
- //assign BCs types and values
- boundaryConditions ( flow );
-
- double center[3];
- for ( int i=0; i<6; i++ ) {
- if ( *flow.boundsIds[i]<0 ) continue;
- CGT::CVector Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
- if ( flow.boundary ( *flow.boundsIds[i] ).useMaxMin ) flow.addBoundingPlane(Normal, *flow.boundsIds[i] );
- else {
- for ( int h=0;h<3;h++ ) center[h] = buffer[*flow.boundsIds[i]].pos[h];
-// cerr << "id="<<*flow.boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl;
- flow.addBoundingPlane ( center, wallThickness, Normal,*flow.boundsIds[i] );
- }
- }
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::triangulate ( Solver& flow )
-{
-///Using Tesselation wrapper (faster)
-// TesselationWrapper TW;
-// if (TW.Tes) delete TW.Tes;
-// TW.Tes = &(flow.tesselation());//point to the current Tes we have in Flowengine
-// TW.insertSceneSpheres();//TW is now really inserting in TemplateFlowEngine, using the faster insert(begin,end)
-// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
-///Using one-by-one insertion
- vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
- FOREACH ( const posData& b, buffer ) {
- if ( !b.exists ) continue;
- if ( b.isSphere ) {
- if (b.id==ignoredBody) continue;
- flow.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
- }
- flow.tesselation().redirected=true;//By inserting one-by-one, we already redirected
- flow.shearLubricationForces.resize ( flow.tesselation().maxId+1 );
- flow.shearLubricationTorques.resize ( flow.tesselation().maxId+1 );
- flow.pumpLubricationTorques.resize ( flow.tesselation().maxId+1 );
- flow.twistLubricationTorques.resize ( flow.tesselation().maxId+1 );
- flow.shearLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
- flow.normalLubricationForce.resize ( flow.tesselation().maxId+1 );
- flow.normalLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::initializeVolumes ( Solver& flow )
-{
- typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
-
- FiniteVerticesIterator vertices_end = flow.tesselation().Triangulation().finite_vertices_end();
- CGT::CVector Zero(0,0,0);
- for (FiniteVerticesIterator V_it = flow.tesselation().Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
-
- FOREACH(CellHandle& cell, flow.tesselation().cellHandles)
- {
- switch ( cell->info().fictious() )
- {
- case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
- case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
- case ( 2 ) : cell->info().volume() = volumeCellDoubleFictious ( cell ); break;
- case ( 3 ) : cell->info().volume() = volumeCellTripleFictious ( cell ); break;
- default: break;
- }
- if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow.volumeSolidPore(cell) ); }
- }
- if (debug) cout << "Volumes initialised." << endl;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageRealCellVelocity()
-{
- solver->averageRelativeCellVelocity();
- Vector3r Vel ( 0,0,0 );
- //AVERAGE CELL VELOCITY
- FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
- for ( int g=0;g<4;g++ ) {
- if ( !cell->vertex ( g )->info().isFictious ) {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- for ( int i=0;i<3;i++ ) Vel[i] = Vel[i] + sph->state->vel[i]/4;
- }
- }
- RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
- CGT::Point pos_av_facet;
- double volume_facet_translation = 0;
- CGT::CVector Vel_av ( Vel[0], Vel[1], Vel[2] );
- for ( int i=0; i<4; i++ ) {
- volume_facet_translation = 0;
- if ( !Tri.is_infinite ( cell->neighbor ( i ) ) ) {
- CGT::CVector Surfk = cell->info()-cell->neighbor ( i )->info();
- Real area = sqrt ( Surfk.squared_length() );
- Surfk = Surfk/area;
- CGT::CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
- pos_av_facet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
- volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
- cell->info().averageVelocity() = cell->info().averageVelocity() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );
- }
- }
- }
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::updateVolumes ( Solver& flow )
-{
- if ( debug ) cout << "Updating volumes.............." << endl;
- Real invDeltaT = 1/scene->dt;
- epsVolMax=0;
- Real totVol=0; Real totDVol=0;
- #ifdef YADE_OPENMP
- const long size=flow.tesselation().cellHandles.size();
- #pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
- for(long i=0; i<size; i++){
- CellHandle& cell = flow.tesselation().cellHandles[i];
- #else
- FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
- #endif
- double newVol, dVol;
- switch ( cell->info().fictious() ) {
- case ( 3 ) : newVol = volumeCellTripleFictious ( cell ); break;
- case ( 2 ) : newVol = volumeCellDoubleFictious ( cell ); break;
- case ( 1 ) : newVol = volumeCellSingleFictious ( cell ); break;
- case ( 0 ) : newVol = volumeCell (cell ); break;
- default: newVol = 0; break;}
- dVol=cell->info().volumeSign* ( newVol - cell->info().volume() );
- cell->info().dv() = dVol*invDeltaT;
- cell->info().volume() = newVol;
- if (defTolerance>0) { //if the criterion is not used, then we skip these updates a save a LOT of time when Nthreads > 1
- #pragma omp atomic
- totVol+=cell->info().volumeSign*newVol;
- #pragma omp atomic
- totDVol+=dVol;}
- }
- if (defTolerance>0) epsVolMax = totDVol/totVol;
- for (unsigned int n=0; n<flow.imposedF.size();n++) {
- flow.IFCells[n]->info().dv()+=flow.imposedF[n].second;
- flow.IFCells[n]->info().Pcondition=false;}
- if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<endl;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellSingleFictious ( Cellhandle cell )
-{
- Vector3r V[3];
- int b=0;
- int w=0;
- cell->info().volumeSign=1;
- Real Wall_coordinate=0;
-
- for ( int y=0;y<4;y++ ) {
- if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
- V[w]=positionBufferCurrent[cell->vertex ( y )->info().id()].pos;
- w++;
- } else {
- b = cell->vertex ( y )->info().id();
- const shared_ptr<Body>& wll = Body::byId ( b , scene );
- if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
- else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
- }
- }
- Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
- return abs ( Volume );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellDoubleFictious ( Cellhandle cell )
-{
- Vector3r A=Vector3r::Zero(), AS=Vector3r::Zero(),B=Vector3r::Zero(), BS=Vector3r::Zero();
-
- cell->info().volumeSign=1;
- int b[2];
- int coord[2];
- Real Wall_coordinate[2];
- int j=0;
- bool first_sph=true;
-
- for ( int g=0;g<4;g++ ) {
- if ( cell->vertex ( g )->info().isFictious ) {
- b[j] = cell->vertex ( g )->info().id();
- coord[j]=solver->boundary ( b[j] ).coordinate;
- if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
- else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
- j++;
- } else if ( first_sph ) {
- A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;
- first_sph=false;
- } else {
- B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;;
- }
- }
- AS[coord[0]]=BS[coord[0]] = Wall_coordinate[0];
-
- //first pyramid with triangular base (A,B,BS)
- Real Vol1=0.5* ( ( A-BS ).cross ( B-BS ) ) [coord[1]]* ( 0.333333333* ( 2*B[coord[1]]+A[coord[1]] )-Wall_coordinate[1] );
- //second pyramid with triangular base (A,AS,BS)
- Real Vol2=0.5* ( ( AS-BS ).cross ( A-BS ) ) [coord[1]]* ( 0.333333333* ( B[coord[1]]+2*A[coord[1]] )-Wall_coordinate[1] );
- return abs ( Vol1+Vol2 );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellTripleFictious ( Cellhandle cell )
-{
- Vector3r A;
-
- int b[3];
- int coord[3];
- Real Wall_coordinate[3];
- int j=0;
- cell->info().volumeSign=1;
-
- for ( int g=0;g<4;g++ ) {
- if ( cell->vertex ( g )->info().isFictious ) {
- b[j] = cell->vertex ( g )->info().id();
- coord[j]=solver->boundary ( b[j] ).coordinate;
- const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
- if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = wll->state->pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
- else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
- j++;
- } else {
- const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
- A= ( sph->state->pos );
- }
- }
- Real Volume = ( A[coord[0]]-Wall_coordinate[0] ) * ( A[coord[1]]-Wall_coordinate[1] ) * ( A[coord[2]]-Wall_coordinate[2] );
- return abs ( Volume );
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-template<class Cellhandle>
-Real TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCell ( Cellhandle cell )
-{
- static const Real inv6 = 1/6.;
- const Vector3r& p0 = positionBufferCurrent[cell->vertex ( 0 )->info().id()].pos;
- const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos;
- const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos;
- const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos;
- Real volume = inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
- if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
- return volume;
-}
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-void TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::computeViscousForces ( Solver& flow )
-{
- if (normalLubrication || shearLubrication || viscousShear){
- if ( debug ) cout << "Application of viscous forces" << endl;
- if ( debug ) cout << "Number of edges = " << flow.edgeIds.size() << endl;
- for ( unsigned int k=0; k<flow.shearLubricationForces.size(); k++ ) flow.shearLubricationForces[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.shearLubricationTorques.size(); k++ ) flow.shearLubricationTorques[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.pumpLubricationTorques.size(); k++ ) flow.pumpLubricationTorques[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.twistLubricationTorques.size(); k++ ) flow.twistLubricationTorques[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.shearLubricationBodyStress.size(); k++) flow.shearLubricationBodyStress[k]=Matrix3r::Zero();
- for ( unsigned int k=0; k<flow.normalLubricationForce.size(); k++ ) flow.normalLubricationForce[k]=Vector3r::Zero();
- for ( unsigned int k=0; k<flow.normalLubricationBodyStress.size(); k++) flow.normalLubricationBodyStress[k]=Matrix3r::Zero();
-
- typedef typename Solver::Tesselation Tesselation;
- const Tesselation& Tes = flow.tesselation();
- flow.deltaShearVel.clear(); flow.normalV.clear(); flow.deltaNormVel.clear(); flow.surfaceDistance.clear(); flow.onlySpheresInteractions.clear(); flow.normalStressInteraction.clear(); flow.shearStressInteraction.clear();
-
-
- for ( int i=0; i< ( int ) flow.edgeIds.size(); i++ ) {
- const VertexInfo& vi1 = *flow.edgeIds[i].first;
- const VertexInfo& vi2 = *flow.edgeIds[i].second;
- const int& id1 = vi1.id();
- const int& id2 = vi2.id();
-
- int hasFictious= Tes.vertex ( id1 )->info().isFictious + Tes.vertex ( id2 )->info().isFictious;
- if (hasFictious>0 or id1==id2) continue;
- const shared_ptr<Body>& sph1 = Body::byId ( id1, scene );
- const shared_ptr<Body>& sph2 = Body::byId ( id2, scene );
- Sphere* s1=YADE_CAST<Sphere*> ( sph1->shape.get() );
- Sphere* s2=YADE_CAST<Sphere*> ( sph2->shape.get() );
- const Real& r1 = s1->radius;
- const Real& r2 = s2->radius;
- Vector3r deltaV; Real deltaNormV; Vector3r deltaShearV;
- Vector3r O1O2Vector; Real O1O2; Vector3r normal; Real surfaceDist; Vector3r O1CVector; Vector3r O2CVector;Real meanRad ;Real Rh; Vector3r deltaAngVel; Vector3r deltaShearAngVel;
- Vector3r shearLubF; Vector3r normaLubF; Vector3r pumpT; Vector3r deltaAngNormVel; Vector3r twistT; Vector3r angVel1; Vector3r angVel2;
- //FIXME: if periodic and velGrad!=0, then deltaV should account for velGrad, not the case currently
- if ( !hasFictious ){
- O1O2Vector = sph2->state->pos + makeVector3r(vi2.ghostShift()) - sph1->state->pos - makeVector3r(vi1.ghostShift());
- O1O2 = O1O2Vector.norm();
- normal= (O1O2Vector/O1O2);
- surfaceDist = O1O2 - r2 - r1;
- O1CVector = (O1O2/2. + (pow(r1,2) - pow(r2,2)) / (2.*O1O2))*normal;
- O2CVector = -(O1O2Vector - O1CVector);
- meanRad = (r2 + r1)/2.;
- Rh = (r1 < r2)? surfaceDist + 0.45 * r1 : surfaceDist + 0.45 * r2;
- deltaV = (sph2->state->vel + sph2->state->angVel.cross(-r2 * normal)) - (sph1->state->vel+ sph1->state->angVel.cross(r1 * normal));
- angVel1 = sph1->state->angVel;
- angVel2 = sph2->state->angVel;
- deltaAngVel = sph2->state->angVel - sph1->state->angVel;
-
- } else {
- if ( hasFictious==1 ) {//for the fictious sphere, use velocity of the boundary, not of the body
- bool v1fictious = Tes.vertex ( id1 )->info().isFictious;
- int bnd = v1fictious? id1 : id2;
- int coord = flow.boundary(bnd).coordinate;
- O1O2 = v1fictious ? abs((sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]) : abs((sph1->state->pos + makeVector3r(Tes.vertex(id1)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]);
- if(v1fictious)
- normal = makeVector3r(flow.boundary(id1).normal);
- else
- normal = -makeVector3r(flow.boundary(id2).normal);
- O1O2Vector = O1O2 * normal;
- meanRad = v1fictious ? r2:r1;
- surfaceDist = O1O2 - meanRad;
- if (v1fictious){
- O1CVector = Vector3r::Zero();
- O2CVector = - O1O2Vector;}
- else{
- O1CVector = O1O2Vector;
- O2CVector = Vector3r::Zero();}
-
- Rh = surfaceDist + 0.45 * meanRad;
- Vector3r v1 = ( Tes.vertex ( id1 )->info().isFictious ) ? flow.boundary ( id1 ).velocity:sph1->state->vel + sph1->state->angVel.cross(r1 * normal);
- Vector3r v2 = ( Tes.vertex ( id2 )->info().isFictious ) ? flow.boundary ( id2 ).velocity:sph2->state->vel + sph2->state->angVel.cross(-r2 * (normal));
- deltaV = v2-v1;
- angVel1 = ( Tes.vertex ( id1 )->info().isFictious ) ? Vector3r::Zero() : sph1->state->angVel;
- angVel2 = ( Tes.vertex ( id2 )->info().isFictious ) ? Vector3r::Zero() : sph2->state->angVel;
- deltaAngVel = angVel2 - angVel1;
- }
- }
- deltaShearV = deltaV - ( normal.dot ( deltaV ) ) *normal;
- deltaShearAngVel = deltaAngVel - ( normal.dot ( deltaAngVel ) ) *normal;
- flow.deltaShearVel.push_back(deltaShearV);
- flow.normalV.push_back(normal);
- flow.surfaceDistance.push_back(max(surfaceDist, 0.) + eps*meanRad);
-
- /// Compute the shear Lubrication force and torque on each particle
-
- if (shearLubrication)
- shearLubF = flow.computeShearLubricationForce(deltaShearV,surfaceDist,i,eps,O1O2,meanRad);
- else if (viscousShear)
- shearLubF = flow.computeViscousShearForce ( deltaShearV, i , Rh);
-
- if (viscousShear || shearLubrication){
-
- flow.shearLubricationForces[id1]+=shearLubF;
- flow.shearLubricationForces[id2]+=(-shearLubF);
- flow.shearLubricationTorques[id1]+=O1CVector.cross(shearLubF);
- flow.shearLubricationTorques[id2]+=O2CVector.cross(-shearLubF);
-
- /// Compute the pump Lubrication torque on each particle
-
- if (pumpTorque){
- pumpT = flow.computePumpTorque(deltaShearAngVel, surfaceDist, i, eps, meanRad );
- flow.pumpLubricationTorques[id1]+=(-pumpT);
- flow.pumpLubricationTorques[id2]+=pumpT;}
-
- /// Compute the twist Lubrication torque on each particle
-
- if (twistTorque){
- deltaAngNormVel = (normal.dot(deltaAngVel))*normal ;
- twistT = flow.computeTwistTorque(deltaAngNormVel, surfaceDist, i, eps, meanRad );
- flow.twistLubricationTorques[id1]+=(-twistT);
- flow.twistLubricationTorques[id2]+=twistT;
- }
- }
- /// Compute the viscous shear stress on each particle
-
- if (viscousShearBodyStress){
- flow.shearLubricationBodyStress[id1] += shearLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
- flow.shearLubricationBodyStress[id2] += (-shearLubF) * O2CVector.transpose()/ (4.0/3.0 *3.14* pow(r2,3));
- flow.shearStressInteraction.push_back(shearLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
- }
-
- /// Compute the normal lubrication force applied on each particle
-
- if (normalLubrication){
- deltaNormV = normal.dot(deltaV);
- flow.deltaNormVel.push_back(deltaNormV * normal);
- normaLubF = flow.computeNormalLubricationForce (deltaNormV, surfaceDist, i,eps,stiffness,scene->dt,meanRad)*normal;
- flow.normalLubricationForce[id1]+=normaLubF;
- flow.normalLubricationForce[id2]+=(-normaLubF);
-
- /// Compute the normal lubrication stress on each particle
-
- if (viscousNormalBodyStress){
- flow.normalLubricationBodyStress[id1] += normaLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
- flow.normalLubricationBodyStress[id2] += (-normaLubF) *O2CVector.transpose() / (4.0/3.0 *3.14* pow(r2,3));
- flow.normalStressInteraction.push_back(normaLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
- }
- }
-
- if (!hasFictious)
- flow.onlySpheresInteractions.push_back(i);
-
- }
- }
-}
-
-template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
-Vector3r TemplateFlowEngine<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageVelocity()
-{
- solver->averageRelativeCellVelocity();
- Vector3r meanVel ( 0,0,0 );
- Real volume=0;
- FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
- for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
- //We could also define velocity using cell's center
-// if ( !cell->info().isReal() ) continue;
- if ( cell->info().isGhost ) continue;
- for ( int i=0;i<3;i++ )
- meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
- volume+=abs ( cell->info().volume() );
- }
- return ( meanVel/volume );
-}
-
-#endif //FLOW_ENGINE
-
-#endif /* YADE_CGAL */
=== added file 'pkg/pfv/FlowEngine.ipp.in'
--- pkg/pfv/FlowEngine.ipp.in 1970-01-01 00:00:00 +0000
+++ pkg/pfv/FlowEngine.ipp.in 2014-06-20 15:28:21 +0000
@@ -0,0 +1,737 @@
+/*************************************************************************
+* Copyright (C) 2009 by Emanuele Catalano <catalano@xxxxxxxxxxxxxxx> *
+* Copyright (C) 2009 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. *
+*************************************************************************/
+#ifdef YADE_CGAL
+
+#ifdef FLOW_ENGINE
+
+#ifdef LINSOLV
+#include <cholmod.h>
+#endif
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::~TemplateFlowEngine_@TEMPLATE_FLOW_NAME@() {}
+
+// YADE_PLUGIN((TFlowEng));
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+unsigned int TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposePressure(Vector3r pos, Real p)
+{
+// if (!flow) LOG_ERROR("no flow defined yet, run at least one iter");
+ solver->imposedP.push_back( pair<CGT::Point,Real>(CGT::Point(pos[0],pos[1],pos[2]),p) );
+ //force immediate update of boundary conditions
+ updateTriangulation=true;
+ return solver->imposedP.size()-1;
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::action()
+{
+ if ( !isActivated ) return;
+ timingDeltas->start();
+ setPositionsBuffer(true);
+ timingDeltas->checkpoint ( "Position buffer" );
+ if (first) {
+ if (multithread) setPositionsBuffer(false);
+ buildTriangulation(pZero,*solver);
+ initializeVolumes(*solver);
+ backgroundSolver=solver;
+ backgroundCompleted=true;
+ }
+ solver->ompThreads = ompThreads>0? ompThreads : omp_get_max_threads();
+
+ timingDeltas->checkpoint ( "Triangulating" );
+ updateVolumes ( *solver );
+ timingDeltas->checkpoint ( "Update_Volumes" );
+
+ epsVolCumulative += epsVolMax;
+ retriangulationLastIter++;
+ if (!updateTriangulation) updateTriangulation = // If not already set true by another function of by the user, check conditions
+ (defTolerance>0 && epsVolCumulative > defTolerance) || retriangulationLastIter>meshUpdateInterval;
+
+ ///compute flow and and forces here
+ if (pressureForce){
+ solver->gaussSeidel(scene->dt);
+ timingDeltas->checkpoint ( "Gauss-Seidel (includes matrix construct and factorization in single-thread mode)" );
+ solver->computeFacetForcesWithCache();}
+ timingDeltas->checkpoint ( "compute_Forces" );
+ ///Application of vicscous forces
+ scene->forces.sync();
+ timingDeltas->checkpoint ( "forces.sync()" );
+ computeViscousForces ( *solver );
+ timingDeltas->checkpoint ( "viscous forces" );
+ Vector3r force;
+ Vector3r torque;
+ FiniteVerticesIterator verticesEnd = solver->T[solver->currentTes].Triangulation().finite_vertices_end();
+ for ( FiniteVerticesIterator vIt = solver->T[solver->currentTes].Triangulation().finite_vertices_begin(); vIt != verticesEnd; vIt++ ) {
+ force = pressureForce ? Vector3r ( vIt->info().forces[0],vIt->info().forces[1],vIt->info().forces[2] ): Vector3r(0,0,0);
+ torque = Vector3r(0,0,0);
+ if (shearLubrication || viscousShear){
+ force = force + solver->shearLubricationForces[vIt->info().id()];
+ torque = torque + solver->shearLubricationTorques[vIt->info().id()];
+ if (pumpTorque)
+ torque = torque + solver->pumpLubricationTorques[vIt->info().id()];
+ }
+ if (twistTorque)
+ torque = torque + solver->twistLubricationTorques[vIt->info().id()];
+ if (normalLubrication)
+ force = force + solver-> normalLubricationForce[vIt->info().id()];
+ scene->forces.addForce ( vIt->info().id(), force);
+ scene->forces.addTorque ( vIt->info().id(), torque);
+ }
+ ///End compute flow and forces
+ timingDeltas->checkpoint ( "Applying Forces" );
+ #ifdef LINSOLV
+ int sleeping = 0;
+ if (multithread && !first) {
+ while (updateTriangulation && !backgroundCompleted) { /*cout<<"sleeping..."<<sleeping++<<endl;*/
+ sleeping++;
+ boost::this_thread::sleep(boost::posix_time::microseconds(1000));}
+ if (debug && sleeping) cerr<<"sleeping..."<<sleeping<<endl;
+ if (updateTriangulation || (ellapsedIter>(0.5*meshUpdateInterval) && backgroundCompleted)) {
+ if (debug) cerr<<"switch flow solver"<<endl;
+ if (useSolver==0) LOG_ERROR("background calculations not available for Gauss-Seidel");
+ if (fluidBulkModulus>0 || doInterpolate) solver->interpolate (solver->T[solver->currentTes], backgroundSolver->T[backgroundSolver->currentTes]);
+ solver=backgroundSolver;
+ backgroundSolver = shared_ptr<FlowSolver> (new FlowSolver);
+ if (metisForced) {backgroundSolver->eSolver.cholmod().nmethods=1; backgroundSolver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;}
+ //Copy imposed pressures/flow from the old solver
+ backgroundSolver->imposedP = vector<pair<CGT::Point,Real> >(solver->imposedP);
+ backgroundSolver->imposedF = vector<pair<CGT::Point,Real> >(solver->imposedF);
+ if (debug) cerr<<"switched"<<endl;
+ setPositionsBuffer(false);//set "parallel" buffer for background calculation
+ backgroundCompleted=false;
+ retriangulationLastIter=ellapsedIter;
+ updateTriangulation=false;
+ epsVolCumulative=0;
+ ellapsedIter=0;
+ boost::thread workerThread(&TemplateFlowEngine_@TEMPLATE_FLOW_NAME@::backgroundAction,this);
+ workerThread.detach();
+ if (debug) cerr<<"backgrounded"<<endl;
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ if (debug) cerr<<"volumes initialized"<<endl;
+ }
+ else {
+ if (debug && !backgroundCompleted) cerr<<"still computing solver in the background, ellapsedIter="<<ellapsedIter<<endl;
+ ellapsedIter++;
+ }
+ } else
+ #endif
+ {
+ if (updateTriangulation && !first) {
+ buildTriangulation (pZero, *solver);
+ initializeVolumes(*solver);
+ computeViscousForces(*solver);
+ updateTriangulation = false;
+ epsVolCumulative=0;
+ retriangulationLastIter=0;
+ ReTrg++;}
+ }
+ first=false;
+ timingDeltas->checkpoint ( "triangulate + init volumes" );
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::backgroundAction()
+{
+ if (useSolver<1) {LOG_ERROR("background calculations not available for Gauss-Seidel"); return;}
+ buildTriangulation ( pZero,*backgroundSolver );
+ //FIXME: GS is computing too much, we need only matrix factorization in fact
+ backgroundSolver->gaussSeidel(scene->dt);
+ //FIXME(2): and here we need only cached variables, not forces
+ backgroundSolver->computeFacetForcesWithCache(/*onlyCache?*/ true);
+// boost::this_thread::sleep(boost::posix_time::seconds(5));
+ backgroundCompleted = true;
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::boundaryConditions ( Solver& flow )
+{
+ for (int k=0;k<6;k++) {
+ flow.boundary (wallIds[k]).flowCondition=!bndCondIsPressure[k];
+ flow.boundary (wallIds[k]).value=bndCondValue[k];
+ flow.boundary (wallIds[k]).velocity = boundaryVelocity[k];//FIXME: needs correct implementation, maybe update the cached pos/vel?
+ }
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::setImposedPressure ( unsigned int cond, Real p)
+{
+ if ( cond>=solver->imposedP.size() ) LOG_ERROR ( "Setting p with cond higher than imposedP size." );
+ solver->imposedP[cond].second=p;
+ //force immediate update of boundary conditions
+ solver->pressureChanged=true;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposeFlux ( Vector3r pos, Real flux){
+ solver->imposedF.push_back ( pair<CGT::Point,Real> ( CGT::Point ( pos[0],pos[1],pos[2] ), flux ) );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedPressure () { solver->imposedP.clear(); solver->IPCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::clearImposedFlux () { solver->imposedF.clear(); solver->IFCells.clear();}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+Real TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::getCellFlux ( unsigned int cond)
+{
+ if ( cond>=solver->imposedP.size() ) {LOG_ERROR ( "Getting flux with cond higher than imposedP size." ); return 0;}
+ double flux=0;
+ typename Solver::CellHandle& cell= solver->IPCells[cond];
+ for ( int ngb=0;ngb<4;ngb++ ) {
+ /*if (!cell->neighbor(ngb)->info().Pcondition)*/
+ flux+= cell->info().kNorm() [ngb]* ( cell->info().p()-cell->neighbor ( ngb )->info().p() );
+ }
+ return flux+cell->info().dv();
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::initSolver ( FlowSolver& flow )
+{
+ flow.Vtotalissimo=0; flow.VSolidTot=0; flow.vPoral=0; flow.sSolidTot=0;
+ flow.slipBoundary=slipBoundary;
+ flow.kFactor = permeabilityFactor;
+ flow.debugOut = debug;
+ flow.useSolver = useSolver;
+ #ifdef EIGENSPARSE_LIB
+ flow.numSolveThreads = numSolveThreads;
+ flow.numFactorizeThreads = numFactorizeThreads;
+ #endif
+ flow.meanKStat = meanKStat;
+ flow.viscosity = viscosity;
+ flow.tolerance=tolerance;
+ flow.relax=relax;
+ flow.clampKValues = clampKValues;
+ flow.maxKdivKmean = maxKdivKmean;
+ flow.minKdivKmean = minKdivKmean;
+ flow.meanKStat = meanKStat;
+ flow.permeabilityMap = permeabilityMap;
+ flow.fluidBulkModulus = fluidBulkModulus;
+// flow.tesselation().Clear();
+ flow.tesselation().maxId=-1;
+ flow.xMin = 1000.0, flow.xMax = -10000.0, flow.yMin = 1000.0, flow.yMax = -10000.0, flow.zMin = 1000.0, flow.zMax = -10000.0;
+}
+
+#ifdef LINSOLV
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::setForceMetis ( bool force )
+{
+ if (force) {
+ metisForced=true;
+ solver->eSolver.cholmod().nmethods=1;
+ solver->eSolver.cholmod().method[0].ordering=CHOLMOD_METIS;
+ } else {cholmod_defaults(&(solver->eSolver.cholmod())); metisForced=false;}
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+bool TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::getForceMetis () {return (solver->eSolver.cholmod().nmethods==1);}
+#endif
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( Solver& flow )
+{
+ buildTriangulation ( 0.f,flow );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::buildTriangulation ( double pZero, Solver& flow )
+{
+ if (first) flow.currentTes=0;
+ else { flow.currentTes=!flow.currentTes; if (debug) cout << "--------RETRIANGULATION-----------" << endl;}
+ flow.resetNetwork();
+ initSolver(flow);
+
+ addBoundary ( flow );
+ triangulate ( flow );
+ if ( debug ) cout << endl << "Tesselating------" << endl << endl;
+ flow.tesselation().compute();
+
+ flow.defineFictiousCells();
+ // For faster loops on cells define this vector
+ flow.tesselation().cellHandles.clear();
+ flow.tesselation().cellHandles.reserve(flow.tesselation().Triangulation().number_of_finite_cells());
+ FiniteCellsIterator cell_end = flow.tesselation().Triangulation().finite_cells_end();
+ int k=0;
+ for ( FiniteCellsIterator cell = flow.tesselation().Triangulation().finite_cells_begin(); cell != cell_end; cell++ ){
+ flow.tesselation().cellHandles.push_back(cell);
+ cell->info().id=k++;}//define unique numbering now, corresponds to position in cellHandles
+ flow.displayStatistics ();
+ flow.computePermeability();
+ //This virtual function does nothing yet, derived class may overload it to make permeability different (see DFN engine)
+ trickPermeability();
+ porosity = flow.vPoralPorosity/flow.vTotalPorosity;
+
+ boundaryConditions ( flow );
+ flow.initializePressure ( pZero );
+
+ if ( !first && !multithread && (useSolver==0 || fluidBulkModulus>0 || doInterpolate)) flow.interpolate ( flow.T[!flow.currentTes], flow.tesselation() );
+ if ( waveAction ) flow.applySinusoidalPressure ( flow.tesselation().Triangulation(), sineMagnitude, sineAverage, 30 );
+ else if (boundaryPressure.size()!=0) flow.applyUserDefinedPressure ( flow.tesselation().Triangulation(), boundaryXPos , boundaryPressure);
+ if (normalLubrication || shearLubrication || viscousShear) flow.computeEdgesSurfaces();
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::setPositionsBuffer(bool current)
+{
+ vector<posData>& buffer = current? positionBufferCurrent : positionBufferParallel;
+ 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;
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::addBoundary ( Solver& flow )
+{
+ vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
+ solver->xMin = Mathr::MAX_REAL, solver->xMax = -Mathr::MAX_REAL, solver->yMin = Mathr::MAX_REAL, solver->yMax = -Mathr::MAX_REAL, solver->zMin = Mathr::MAX_REAL, solver->zMax = -Mathr::MAX_REAL;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists ) continue;
+ if ( b.isSphere ) {
+ const Real& rad = b.radius;
+ const Real& x = b.pos[0];
+ const Real& y = b.pos[1];
+ const Real& z = b.pos[2];
+ flow.xMin = min ( flow.xMin, x-rad );
+ flow.xMax = max ( flow.xMax, x+rad );
+ flow.yMin = min ( flow.yMin, y-rad );
+ flow.yMax = max ( flow.yMax, y+rad );
+ flow.zMin = min ( flow.zMin, z-rad );
+ flow.zMax = max ( flow.zMax, z+rad );
+ }
+ }
+ //FIXME idOffset must be set correctly, not the case here (always 0), then we need walls first or it will fail
+ idOffset = flow.tesselation().maxId+1;
+ flow.idOffset = idOffset;
+ flow.sectionArea = ( flow.xMax - flow.xMin ) * ( flow.zMax-flow.zMin );
+ flow.vTotal = ( flow.xMax-flow.xMin ) * ( flow.yMax-flow.yMin ) * ( flow.zMax-flow.zMin );
+ flow.yMinId=wallIds[ymin];
+ flow.yMaxId=wallIds[ymax];
+ flow.xMaxId=wallIds[xmax];
+ flow.xMinId=wallIds[xmin];
+ flow.zMinId=wallIds[zmin];
+ flow.zMaxId=wallIds[zmax];
+
+ //FIXME: Id's order in boundsIds is done according to the enumeration of boundaries from TXStressController.hpp, line 31. DON'T CHANGE IT!
+ flow.boundsIds[0]= &flow.xMinId;
+ flow.boundsIds[1]= &flow.xMaxId;
+ flow.boundsIds[2]= &flow.yMinId;
+ flow.boundsIds[3]= &flow.yMaxId;
+ flow.boundsIds[4]= &flow.zMinId;
+ flow.boundsIds[5]= &flow.zMaxId;
+
+ for (int k=0;k<6;k++) flow.boundary ( *flow.boundsIds[k] ).useMaxMin = boundaryUseMaxMin[k];
+
+ flow.cornerMin = CGT::Point ( flow.xMin, flow.yMin, flow.zMin );
+ flow.cornerMax = CGT::Point ( flow.xMax, flow.yMax, flow.zMax );
+
+ //assign BCs types and values
+ boundaryConditions ( flow );
+
+ double center[3];
+ for ( int i=0; i<6; i++ ) {
+ if ( *flow.boundsIds[i]<0 ) continue;
+ CGT::CVector Normal ( normal[i].x(), normal[i].y(), normal[i].z() );
+ if ( flow.boundary ( *flow.boundsIds[i] ).useMaxMin ) flow.addBoundingPlane(Normal, *flow.boundsIds[i] );
+ else {
+ for ( int h=0;h<3;h++ ) center[h] = buffer[*flow.boundsIds[i]].pos[h];
+// cerr << "id="<<*flow.boundsIds[i] <<" center="<<center[0]<<","<<center[1]<<","<<center[2]<<endl;
+ flow.addBoundingPlane ( center, wallThickness, Normal,*flow.boundsIds[i] );
+ }
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::triangulate ( Solver& flow )
+{
+///Using Tesselation wrapper (faster)
+// TesselationWrapper TW;
+// if (TW.Tes) delete TW.Tes;
+// TW.Tes = &(flow.tesselation());//point to the current Tes we have in Flowengine
+// TW.insertSceneSpheres();//TW is now really inserting in TemplateFlowEngine_@TEMPLATE_FLOW_NAME@, using the faster insert(begin,end)
+// TW.Tes = NULL;//otherwise, Tes would be deleted by ~TesselationWrapper() at the end of the function.
+///Using one-by-one insertion
+ vector<posData>& buffer = multithread ? positionBufferParallel : positionBufferCurrent;
+ FOREACH ( const posData& b, buffer ) {
+ if ( !b.exists ) continue;
+ if ( b.isSphere ) {
+ if (b.id==ignoredBody) continue;
+ flow.tesselation().insert ( b.pos[0], b.pos[1], b.pos[2], b.radius, b.id );}
+ }
+ flow.tesselation().redirected=true;//By inserting one-by-one, we already redirected
+ flow.shearLubricationForces.resize ( flow.tesselation().maxId+1 );
+ flow.shearLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.pumpLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.twistLubricationTorques.resize ( flow.tesselation().maxId+1 );
+ flow.shearLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
+ flow.normalLubricationForce.resize ( flow.tesselation().maxId+1 );
+ flow.normalLubricationBodyStress.resize ( flow.tesselation().maxId+1 );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::initializeVolumes ( Solver& flow )
+{
+ typedef typename Solver::FiniteVerticesIterator FiniteVerticesIterator;
+
+ FiniteVerticesIterator vertices_end = flow.tesselation().Triangulation().finite_vertices_end();
+ CGT::CVector Zero(0,0,0);
+ for (FiniteVerticesIterator V_it = flow.tesselation().Triangulation().finite_vertices_begin(); V_it!= vertices_end; V_it++) V_it->info().forces=Zero;
+
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles)
+ {
+ switch ( cell->info().fictious() )
+ {
+ case ( 0 ) : cell->info().volume() = volumeCell ( cell ); break;
+ case ( 1 ) : cell->info().volume() = volumeCellSingleFictious ( cell ); break;
+ case ( 2 ) : cell->info().volume() = volumeCellDoubleFictious ( cell ); break;
+ case ( 3 ) : cell->info().volume() = volumeCellTripleFictious ( cell ); break;
+ default: break;
+ }
+ if (flow.fluidBulkModulus>0) { cell->info().invVoidVolume() = 1 / ( abs(cell->info().volume()) - flow.volumeSolidPore(cell) ); }
+ }
+ if (debug) cout << "Volumes initialised." << endl;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageRealCellVelocity()
+{
+ solver->averageRelativeCellVelocity();
+ Vector3r Vel ( 0,0,0 );
+ //AVERAGE CELL VELOCITY
+ FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
+ for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ for ( int g=0;g<4;g++ ) {
+ if ( !cell->vertex ( g )->info().isFictious ) {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ for ( int i=0;i<3;i++ ) Vel[i] = Vel[i] + sph->state->vel[i]/4;
+ }
+ }
+ RTriangulation& Tri = solver->T[solver->currentTes].Triangulation();
+ CGT::Point pos_av_facet;
+ double volume_facet_translation = 0;
+ CGT::CVector Vel_av ( Vel[0], Vel[1], Vel[2] );
+ for ( int i=0; i<4; i++ ) {
+ volume_facet_translation = 0;
+ if ( !Tri.is_infinite ( cell->neighbor ( i ) ) ) {
+ CGT::CVector Surfk = cell->info()-cell->neighbor ( i )->info();
+ Real area = sqrt ( Surfk.squared_length() );
+ Surfk = Surfk/area;
+ CGT::CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+ pos_av_facet = ( CGT::Point ) cell->info() + ( branch*Surfk ) *Surfk;
+ volume_facet_translation += Vel_av*cell->info().facetSurfaces[i];
+ cell->info().averageVelocity() = cell->info().averageVelocity() - volume_facet_translation/cell->info().volume() * ( pos_av_facet-CGAL::ORIGIN );
+ }
+ }
+ }
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::updateVolumes ( Solver& flow )
+{
+ if ( debug ) cout << "Updating volumes.............." << endl;
+ Real invDeltaT = 1/scene->dt;
+ epsVolMax=0;
+ Real totVol=0; Real totDVol=0;
+ #ifdef YADE_OPENMP
+ const long size=flow.tesselation().cellHandles.size();
+ #pragma omp parallel for num_threads(ompThreads>0 ? ompThreads : 1)
+ for(long i=0; i<size; i++){
+ CellHandle& cell = flow.tesselation().cellHandles[i];
+ #else
+ FOREACH(CellHandle& cell, flow.tesselation().cellHandles){
+ #endif
+ double newVol, dVol;
+ switch ( cell->info().fictious() ) {
+ case ( 3 ) : newVol = volumeCellTripleFictious ( cell ); break;
+ case ( 2 ) : newVol = volumeCellDoubleFictious ( cell ); break;
+ case ( 1 ) : newVol = volumeCellSingleFictious ( cell ); break;
+ case ( 0 ) : newVol = volumeCell (cell ); break;
+ default: newVol = 0; break;}
+ dVol=cell->info().volumeSign* ( newVol - cell->info().volume() );
+ cell->info().dv() = dVol*invDeltaT;
+ cell->info().volume() = newVol;
+ if (defTolerance>0) { //if the criterion is not used, then we skip these updates a save a LOT of time when Nthreads > 1
+ #pragma omp atomic
+ totVol+=cell->info().volumeSign*newVol;
+ #pragma omp atomic
+ totDVol+=dVol;}
+ }
+ if (defTolerance>0) epsVolMax = totDVol/totVol;
+ for (unsigned int n=0; n<flow.imposedF.size();n++) {
+ flow.IFCells[n]->info().dv()+=flow.imposedF[n].second;
+ flow.IFCells[n]->info().Pcondition=false;}
+ if ( debug ) cout << "Updated volumes, total =" <<totVol<<", dVol="<<totDVol<<endl;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellSingleFictious ( Cellhandle cell )
+{
+ Vector3r V[3];
+ int b=0;
+ int w=0;
+ cell->info().volumeSign=1;
+ Real Wall_coordinate=0;
+
+ for ( int y=0;y<4;y++ ) {
+ if ( ! ( cell->vertex ( y )->info().isFictious ) ) {
+ V[w]=positionBufferCurrent[cell->vertex ( y )->info().id()].pos;
+ w++;
+ } else {
+ b = cell->vertex ( y )->info().id();
+ const shared_ptr<Body>& wll = Body::byId ( b , scene );
+ if ( !solver->boundary ( b ).useMaxMin ) Wall_coordinate = wll->state->pos[solver->boundary ( b ).coordinate]+ ( solver->boundary ( b ).normal[solver->boundary ( b ).coordinate] ) *wallThickness/2.;
+ else Wall_coordinate = solver->boundary ( b ).p[solver->boundary ( b ).coordinate];
+ }
+ }
+ Real Volume = 0.5* ( ( V[0]-V[1] ).cross ( V[0]-V[2] ) ) [solver->boundary ( b ).coordinate] * ( 0.33333333333* ( V[0][solver->boundary ( b ).coordinate]+ V[1][solver->boundary ( b ).coordinate]+ V[2][solver->boundary ( b ).coordinate] ) - Wall_coordinate );
+ return abs ( Volume );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellDoubleFictious ( Cellhandle cell )
+{
+ Vector3r A=Vector3r::Zero(), AS=Vector3r::Zero(),B=Vector3r::Zero(), BS=Vector3r::Zero();
+
+ cell->info().volumeSign=1;
+ int b[2];
+ int coord[2];
+ Real Wall_coordinate[2];
+ int j=0;
+ bool first_sph=true;
+
+ for ( int g=0;g<4;g++ ) {
+ if ( cell->vertex ( g )->info().isFictious ) {
+ b[j] = cell->vertex ( g )->info().id();
+ coord[j]=solver->boundary ( b[j] ).coordinate;
+ if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = positionBufferCurrent[b[j]].pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
+ else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
+ j++;
+ } else if ( first_sph ) {
+ A=AS=/*AT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;
+ first_sph=false;
+ } else {
+ B=BS=/*BT=*/ positionBufferCurrent[cell->vertex(g)->info().id()].pos;;
+ }
+ }
+ AS[coord[0]]=BS[coord[0]] = Wall_coordinate[0];
+
+ //first pyramid with triangular base (A,B,BS)
+ Real Vol1=0.5* ( ( A-BS ).cross ( B-BS ) ) [coord[1]]* ( 0.333333333* ( 2*B[coord[1]]+A[coord[1]] )-Wall_coordinate[1] );
+ //second pyramid with triangular base (A,AS,BS)
+ Real Vol2=0.5* ( ( AS-BS ).cross ( A-BS ) ) [coord[1]]* ( 0.333333333* ( B[coord[1]]+2*A[coord[1]] )-Wall_coordinate[1] );
+ return abs ( Vol1+Vol2 );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCellTripleFictious ( Cellhandle cell )
+{
+ Vector3r A;
+
+ int b[3];
+ int coord[3];
+ Real Wall_coordinate[3];
+ int j=0;
+ cell->info().volumeSign=1;
+
+ for ( int g=0;g<4;g++ ) {
+ if ( cell->vertex ( g )->info().isFictious ) {
+ b[j] = cell->vertex ( g )->info().id();
+ coord[j]=solver->boundary ( b[j] ).coordinate;
+ const shared_ptr<Body>& wll = Body::byId ( b[j] , scene );
+ if ( !solver->boundary ( b[j] ).useMaxMin ) Wall_coordinate[j] = wll->state->pos[coord[j]] + ( solver->boundary ( b[j] ).normal[coord[j]] ) *wallThickness/2.;
+ else Wall_coordinate[j] = solver->boundary ( b[j] ).p[coord[j]];
+ j++;
+ } else {
+ const shared_ptr<Body>& sph = Body::byId ( cell->vertex ( g )->info().id(), scene );
+ A= ( sph->state->pos );
+ }
+ }
+ Real Volume = ( A[coord[0]]-Wall_coordinate[0] ) * ( A[coord[1]]-Wall_coordinate[1] ) * ( A[coord[2]]-Wall_coordinate[2] );
+ return abs ( Volume );
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+template<class Cellhandle>
+Real TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::volumeCell ( Cellhandle cell )
+{
+ static const Real inv6 = 1/6.;
+ const Vector3r& p0 = positionBufferCurrent[cell->vertex ( 0 )->info().id()].pos;
+ const Vector3r& p1 = positionBufferCurrent[cell->vertex ( 1 )->info().id()].pos;
+ const Vector3r& p2 = positionBufferCurrent[cell->vertex ( 2 )->info().id()].pos;
+ const Vector3r& p3 = positionBufferCurrent[cell->vertex ( 3 )->info().id()].pos;
+ Real volume = inv6 * ((p0-p1).cross(p0-p2)).dot(p0-p3);
+ if ( ! ( cell->info().volumeSign ) ) cell->info().volumeSign= ( volume>0 ) ?1:-1;
+ return volume;
+}
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::computeViscousForces ( Solver& flow )
+{
+ if (normalLubrication || shearLubrication || viscousShear){
+ if ( debug ) cout << "Application of viscous forces" << endl;
+ if ( debug ) cout << "Number of edges = " << flow.edgeIds.size() << endl;
+ for ( unsigned int k=0; k<flow.shearLubricationForces.size(); k++ ) flow.shearLubricationForces[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.shearLubricationTorques.size(); k++ ) flow.shearLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.pumpLubricationTorques.size(); k++ ) flow.pumpLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.twistLubricationTorques.size(); k++ ) flow.twistLubricationTorques[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.shearLubricationBodyStress.size(); k++) flow.shearLubricationBodyStress[k]=Matrix3r::Zero();
+ for ( unsigned int k=0; k<flow.normalLubricationForce.size(); k++ ) flow.normalLubricationForce[k]=Vector3r::Zero();
+ for ( unsigned int k=0; k<flow.normalLubricationBodyStress.size(); k++) flow.normalLubricationBodyStress[k]=Matrix3r::Zero();
+
+ typedef typename Solver::Tesselation Tesselation;
+ const Tesselation& Tes = flow.tesselation();
+ flow.deltaShearVel.clear(); flow.normalV.clear(); flow.deltaNormVel.clear(); flow.surfaceDistance.clear(); flow.onlySpheresInteractions.clear(); flow.normalStressInteraction.clear(); flow.shearStressInteraction.clear();
+
+
+ for ( int i=0; i< ( int ) flow.edgeIds.size(); i++ ) {
+ const VertexInfo& vi1 = *flow.edgeIds[i].first;
+ const VertexInfo& vi2 = *flow.edgeIds[i].second;
+ const int& id1 = vi1.id();
+ const int& id2 = vi2.id();
+
+ int hasFictious= Tes.vertex ( id1 )->info().isFictious + Tes.vertex ( id2 )->info().isFictious;
+ if (hasFictious>0 or id1==id2) continue;
+ const shared_ptr<Body>& sph1 = Body::byId ( id1, scene );
+ const shared_ptr<Body>& sph2 = Body::byId ( id2, scene );
+ Sphere* s1=YADE_CAST<Sphere*> ( sph1->shape.get() );
+ Sphere* s2=YADE_CAST<Sphere*> ( sph2->shape.get() );
+ const Real& r1 = s1->radius;
+ const Real& r2 = s2->radius;
+ Vector3r deltaV; Real deltaNormV; Vector3r deltaShearV;
+ Vector3r O1O2Vector; Real O1O2; Vector3r normal; Real surfaceDist; Vector3r O1CVector; Vector3r O2CVector;Real meanRad ;Real Rh; Vector3r deltaAngVel; Vector3r deltaShearAngVel;
+ Vector3r shearLubF; Vector3r normaLubF; Vector3r pumpT; Vector3r deltaAngNormVel; Vector3r twistT; Vector3r angVel1; Vector3r angVel2;
+ //FIXME: if periodic and velGrad!=0, then deltaV should account for velGrad, not the case currently
+ if ( !hasFictious ){
+ O1O2Vector = sph2->state->pos + makeVector3r(vi2.ghostShift()) - sph1->state->pos - makeVector3r(vi1.ghostShift());
+ O1O2 = O1O2Vector.norm();
+ normal= (O1O2Vector/O1O2);
+ surfaceDist = O1O2 - r2 - r1;
+ O1CVector = (O1O2/2. + (pow(r1,2) - pow(r2,2)) / (2.*O1O2))*normal;
+ O2CVector = -(O1O2Vector - O1CVector);
+ meanRad = (r2 + r1)/2.;
+ Rh = (r1 < r2)? surfaceDist + 0.45 * r1 : surfaceDist + 0.45 * r2;
+ deltaV = (sph2->state->vel + sph2->state->angVel.cross(-r2 * normal)) - (sph1->state->vel+ sph1->state->angVel.cross(r1 * normal));
+ angVel1 = sph1->state->angVel;
+ angVel2 = sph2->state->angVel;
+ deltaAngVel = sph2->state->angVel - sph1->state->angVel;
+
+ } else {
+ if ( hasFictious==1 ) {//for the fictious sphere, use velocity of the boundary, not of the body
+ bool v1fictious = Tes.vertex ( id1 )->info().isFictious;
+ int bnd = v1fictious? id1 : id2;
+ int coord = flow.boundary(bnd).coordinate;
+ O1O2 = v1fictious ? abs((sph2->state->pos + makeVector3r(Tes.vertex(id2)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]) : abs((sph1->state->pos + makeVector3r(Tes.vertex(id1)->info().ghostShift()))[coord] - flow.boundary(bnd).p[coord]);
+ if(v1fictious)
+ normal = makeVector3r(flow.boundary(id1).normal);
+ else
+ normal = -makeVector3r(flow.boundary(id2).normal);
+ O1O2Vector = O1O2 * normal;
+ meanRad = v1fictious ? r2:r1;
+ surfaceDist = O1O2 - meanRad;
+ if (v1fictious){
+ O1CVector = Vector3r::Zero();
+ O2CVector = - O1O2Vector;}
+ else{
+ O1CVector = O1O2Vector;
+ O2CVector = Vector3r::Zero();}
+
+ Rh = surfaceDist + 0.45 * meanRad;
+ Vector3r v1 = ( Tes.vertex ( id1 )->info().isFictious ) ? flow.boundary ( id1 ).velocity:sph1->state->vel + sph1->state->angVel.cross(r1 * normal);
+ Vector3r v2 = ( Tes.vertex ( id2 )->info().isFictious ) ? flow.boundary ( id2 ).velocity:sph2->state->vel + sph2->state->angVel.cross(-r2 * (normal));
+ deltaV = v2-v1;
+ angVel1 = ( Tes.vertex ( id1 )->info().isFictious ) ? Vector3r::Zero() : sph1->state->angVel;
+ angVel2 = ( Tes.vertex ( id2 )->info().isFictious ) ? Vector3r::Zero() : sph2->state->angVel;
+ deltaAngVel = angVel2 - angVel1;
+ }
+ }
+ deltaShearV = deltaV - ( normal.dot ( deltaV ) ) *normal;
+ deltaShearAngVel = deltaAngVel - ( normal.dot ( deltaAngVel ) ) *normal;
+ flow.deltaShearVel.push_back(deltaShearV);
+ flow.normalV.push_back(normal);
+ flow.surfaceDistance.push_back(max(surfaceDist, 0.) + eps*meanRad);
+
+ /// Compute the shear Lubrication force and torque on each particle
+
+ if (shearLubrication)
+ shearLubF = flow.computeShearLubricationForce(deltaShearV,surfaceDist,i,eps,O1O2,meanRad);
+ else if (viscousShear)
+ shearLubF = flow.computeViscousShearForce ( deltaShearV, i , Rh);
+
+ if (viscousShear || shearLubrication){
+
+ flow.shearLubricationForces[id1]+=shearLubF;
+ flow.shearLubricationForces[id2]+=(-shearLubF);
+ flow.shearLubricationTorques[id1]+=O1CVector.cross(shearLubF);
+ flow.shearLubricationTorques[id2]+=O2CVector.cross(-shearLubF);
+
+ /// Compute the pump Lubrication torque on each particle
+
+ if (pumpTorque){
+ pumpT = flow.computePumpTorque(deltaShearAngVel, surfaceDist, i, eps, meanRad );
+ flow.pumpLubricationTorques[id1]+=(-pumpT);
+ flow.pumpLubricationTorques[id2]+=pumpT;}
+
+ /// Compute the twist Lubrication torque on each particle
+
+ if (twistTorque){
+ deltaAngNormVel = (normal.dot(deltaAngVel))*normal ;
+ twistT = flow.computeTwistTorque(deltaAngNormVel, surfaceDist, i, eps, meanRad );
+ flow.twistLubricationTorques[id1]+=(-twistT);
+ flow.twistLubricationTorques[id2]+=twistT;
+ }
+ }
+ /// Compute the viscous shear stress on each particle
+
+ if (viscousShearBodyStress){
+ flow.shearLubricationBodyStress[id1] += shearLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
+ flow.shearLubricationBodyStress[id2] += (-shearLubF) * O2CVector.transpose()/ (4.0/3.0 *3.14* pow(r2,3));
+ flow.shearStressInteraction.push_back(shearLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
+ }
+
+ /// Compute the normal lubrication force applied on each particle
+
+ if (normalLubrication){
+ deltaNormV = normal.dot(deltaV);
+ flow.deltaNormVel.push_back(deltaNormV * normal);
+ normaLubF = flow.computeNormalLubricationForce (deltaNormV, surfaceDist, i,eps,stiffness,scene->dt,meanRad)*normal;
+ flow.normalLubricationForce[id1]+=normaLubF;
+ flow.normalLubricationForce[id2]+=(-normaLubF);
+
+ /// Compute the normal lubrication stress on each particle
+
+ if (viscousNormalBodyStress){
+ flow.normalLubricationBodyStress[id1] += normaLubF * O1CVector.transpose()/ (4.0/3.0 *3.14* pow(r1,3));
+ flow.normalLubricationBodyStress[id2] += (-normaLubF) *O2CVector.transpose() / (4.0/3.0 *3.14* pow(r2,3));
+ flow.normalStressInteraction.push_back(normaLubF * O1O2Vector.transpose()/(4.0/3.0 *3.14* pow(r1,3)));
+ }
+ }
+
+ if (!hasFictious)
+ flow.onlySpheresInteractions.push_back(i);
+
+ }
+ }
+}
+
+template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
+Vector3r TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::averageVelocity()
+{
+ solver->averageRelativeCellVelocity();
+ Vector3r meanVel ( 0,0,0 );
+ Real volume=0;
+ FiniteCellsIterator cell_end = solver->T[solver->currentTes].Triangulation().finite_cells_end();
+ for ( FiniteCellsIterator cell = solver->T[solver->currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++ ) {
+ //We could also define velocity using cell's center
+// if ( !cell->info().isReal() ) continue;
+ if ( cell->info().isGhost ) continue;
+ for ( int i=0;i<3;i++ )
+ meanVel[i]=meanVel[i]+ ( ( cell->info().averageVelocity() ) [i] * abs ( cell->info().volume() ) );
+ volume+=abs ( cell->info().volume() );
+ }
+ return ( meanVel/volume );
+}
+
+#endif //FLOW_ENGINE
+
+#endif /* YADE_CGAL */
+
=== modified file 'pkg/pfv/PeriodicFlowEngine.cpp'
--- pkg/pfv/PeriodicFlowEngine.cpp 2014-06-04 17:19:50 +0000
+++ pkg/pfv/PeriodicFlowEngine.cpp 2014-06-20 17:48:57 +0000
@@ -13,10 +13,9 @@
/// It is a bit more complicated as for FlowEngine, though, because we need template inheriting from template, which breaks YADE_CLASS_XXX logic_error
/// See below the commented exemple, for a possible solution
-#define TEMPLATE_FLOW_NAME FlowEngine_PeriodicInfo
-#include <yade/pkg/pfv/FlowEngine.hpp>
+#include "FlowEngine_FlowEngine_PeriodicInfo.hpp"
-class PeriodicCellInfo : public FlowCellInfo
+class PeriodicCellInfo : public FlowCellInfo_FlowEngine_PeriodicInfo
{
public:
static CVector gradP;
@@ -43,7 +42,7 @@
};
-class PeriodicVertexInfo : public FlowVertexInfo {
+class PeriodicVertexInfo : public FlowVertexInfo_FlowEngine_PeriodicInfo {
public:
PeriodicVertexInfo& operator= (const CVector &u) { CVector::operator= (u); return *this; }
PeriodicVertexInfo& operator= (const float &scalar) { s=scalar; return *this; }
@@ -64,7 +63,7 @@
#define _PeriFlowSolver CGT::PeriodicFlow<PeriFlowTesselation>
#endif
-typedef TemplateFlowEngine< PeriodicCellInfo,
+typedef TemplateFlowEngine_FlowEngine_PeriodicInfo< PeriodicCellInfo,
PeriodicVertexInfo,
CGT::PeriodicTesselation<CGT::_Tesselation<CGT::TriangulationTypes<PeriodicVertexInfo,PeriodicCellInfo> > >,
_PeriFlowSolver>
@@ -529,4 +528,4 @@
YADE_PLUGIN((PeriodicFlowEngine));
#endif //FLOW_ENGINE
-#endif /* YADE_CGAL */
\ No newline at end of file
+#endif /* YADE_CGAL */
=== modified file 'pkg/pfv/SoluteFlowEngine.cpp'
--- pkg/pfv/SoluteFlowEngine.cpp 2014-05-23 13:20:43 +0000
+++ pkg/pfv/SoluteFlowEngine.cpp 2014-06-20 20:10:04 +0000
@@ -11,22 +11,21 @@
// #define SOLUTE_FLOW
#ifdef SOLUTE_FLOW
-#define TEMPLATE_FLOW_NAME SoluteFlowEngineT
-#include <yade/pkg/pfv/FlowEngine.hpp>
+#include "FlowEngine_SoluteFlowEngineT.hpp"
#include <Eigen/Sparse>
-class SoluteCellInfo : public FlowCellInfo
+class SoluteCellInfo : public FlowCellInfo_SoluteFlowEngineT
{
public:
Real solute_concentration;
- SoluteCellInfo (void) : FlowCellInfo() {solute_concentration=0;}
+ SoluteCellInfo (void) : FlowCellInfo_SoluteFlowEngineT() {solute_concentration=0;}
inline Real& solute (void) {return solute_concentration;}
inline const Real& solute (void) const {return solute_concentration;}
- inline void getInfo (const SoluteCellInfo& otherCellInfo) {FlowCellInfo::getInfo(otherCellInfo); solute()=otherCellInfo.solute();}
+ inline void getInfo (const SoluteCellInfo& otherCellInfo) {FlowCellInfo_SoluteFlowEngineT::getInfo(otherCellInfo); solute()=otherCellInfo.solute();}
};
-typedef TemplateFlowEngine<SoluteCellInfo,FlowVertexInfo> SoluteFlowEngineT;
+typedef TemplateFlowEngine_SoluteFlowEngineT<SoluteCellInfo,FlowVertexInfo_SoluteFlowEngineT> SoluteFlowEngineT;
REGISTER_SERIALIZABLE(SoluteFlowEngineT);
YADE_PLUGIN((SoluteFlowEngineT));
=== modified file 'py/bodiesHandling.py'
--- py/bodiesHandling.py 2013-03-08 21:26:47 +0000
+++ py/bodiesHandling.py 2014-06-12 16:54:31 +0000
@@ -168,7 +168,7 @@
return dimensions
#spheresPackDimensions==================================================
-def spheresModify(idSpheres=[],mask=-1,shift=Vector3.Zero,scale=1.0,orientation=Quaternion.Identity,copy=False):
+def spheresModify(idSpheres=[],mask=-1,shift=Vector3.Zero,scale=1.0,orientation=Quaternion((0,1,0),0.0),copy=False):
"""The function accepts the list of spheres id's or list of bodies and modifies them: rotating, scaling, shifting.
if copy=True copies bodies and modifies them.
Also the mask can be given. If idSpheres not empty, the function affects only bodies, where the mask passes.
=== modified file 'py/geom.py'
--- py/geom.py 2014-01-14 18:37:09 +0000
+++ py/geom.py 2014-06-12 16:54:31 +0000
@@ -12,7 +12,7 @@
from miniEigen import *
#facetBox===============================================================
-def facetBox(center,extents,orientation=Quaternion.Identity,wallMask=63,**kw):
+def facetBox(center,extents,orientation=Quaternion((0,1,0),0.0),wallMask=63,**kw):
"""
Create arbitrarily-aligned box composed of facets, with given center, extents and orientation.
If any of the box dimensions is zero, corresponding facets will not be created. The facets are oriented outwards from the box.
@@ -28,7 +28,7 @@
return facetParallelepiped(center=center, extents=extents, height=extents[2], orientation=orientation, wallMask=wallMask, **kw)
#facetParallelepiped===============================================================
-def facetParallelepiped(center,extents,height,orientation=Quaternion.Identity,wallMask=63,**kw):
+def facetParallelepiped(center,extents,height,orientation=Quaternion((0,1,0),0.0),wallMask=63,**kw):
"""
Create arbitrarily-aligned Parallelepiped composed of facets, with given center, extents, height and orientation.
If any of the parallelepiped dimensions is zero, corresponding facets will not be created. The facets are oriented outwards from the parallelepiped.
@@ -90,7 +90,7 @@
return ret
#facetCylinder==========================================================
-def facetCylinder(center,radius,height,orientation=Quaternion.Identity,
+def facetCylinder(center,radius,height,orientation=Quaternion((0,1,0),0.0),
segmentsNumber=10,wallMask=7,angleRange=None,closeGap=False,
radiusTopInner=-1, radiusBottomInner=-1,
**kw):
@@ -164,7 +164,7 @@
#facetCone==============================================================
-def facetCone(center,radiusTop,radiusBottom,height,orientation=Quaternion.Identity,
+def facetCone(center,radiusTop,radiusBottom,height,orientation=Quaternion((0,1,0),0.0),
segmentsNumber=10,wallMask=7,angleRange=None,closeGap=False,
radiusTopInner=-1, radiusBottomInner=-1,
**kw):
@@ -195,7 +195,7 @@
**kw)
#facetPolygon===========================================================
-def facetPolygon(center,radiusOuter,orientation=Quaternion.Identity,segmentsNumber=10,angleRange=None,radiusInner=0,**kw):
+def facetPolygon(center,radiusOuter,orientation=Quaternion((0,1,0),0.0),segmentsNumber=10,angleRange=None,radiusInner=0,**kw):
"""
Create arbitrarily-aligned polygon composed of facets, with given center, radius (outer and inner) and orientation.
Return List of facets forming the polygon;
@@ -214,7 +214,7 @@
return facetPolygonHelixGenerator(center=center,radiusOuter=radiusOuter,orientation=orientation,segmentsNumber=segmentsNumber,angleRange=angleRange,radiusInner=radiusInner,**kw)
#facetHelix===========================================================
-def facetHelix(center,radiusOuter,pitch,orientation=Quaternion.Identity,segmentsNumber=10,angleRange=None,radiusInner=0,**kw):
+def facetHelix(center,radiusOuter,pitch,orientation=Quaternion((0,1,0),0.0),segmentsNumber=10,angleRange=None,radiusInner=0,**kw):
"""
Create arbitrarily-aligned helix composed of facets, with given center, radius (outer and inner), pitch and orientation.
Return List of facets forming the helix;
@@ -233,7 +233,7 @@
return facetPolygonHelixGenerator(center=center,radiusOuter=radiusOuter,orientation=orientation,segmentsNumber=segmentsNumber,angleRange=angleRange,radiusInner=radiusInner,pitch=pitch,**kw)
#facetBunker============================================================
-def facetBunker(center,dBunker,dOutput,hBunker,hOutput,hPipe=0.0,orientation=Quaternion.Identity,segmentsNumber=10,wallMask=4,angleRange=None,closeGap=False,**kw):
+def facetBunker(center,dBunker,dOutput,hBunker,hOutput,hPipe=0.0,orientation=Quaternion((0,1,0),0.0),segmentsNumber=10,wallMask=4,angleRange=None,closeGap=False,**kw):
"""
Create arbitrarily-aligned bunker, composed of facets, with given center, radii, heights and orientation.
Return List of facets forming the bunker;
@@ -297,7 +297,7 @@
return ret
#facetPolygonHelixGenerator==================================================
-def facetPolygonHelixGenerator(center,radiusOuter,pitch=0,orientation=Quaternion.Identity,segmentsNumber=10,angleRange=None,radiusInner=0,**kw):
+def facetPolygonHelixGenerator(center,radiusOuter,pitch=0,orientation=Quaternion((0,1,0),0.0),segmentsNumber=10,angleRange=None,radiusInner=0,**kw):
"""
Please, do not use this function directly! Use geom.facetPloygon and geom.facetHelix instead.
This is the base function for generating polygons and helixes from facets.
@@ -342,7 +342,7 @@
#facetCylinderConeGenerator=============================================
-def facetCylinderConeGenerator(center,radiusTop,height,orientation=Quaternion.Identity,
+def facetCylinderConeGenerator(center,radiusTop,height,orientation=Quaternion((0,1,0),0.0),
segmentsNumber=10,wallMask=7,angleRange=None,closeGap=False,
radiusBottom=-1,
radiusTopInner=-1,
=== modified file 'py/wrapper/customConverters.cpp'
--- py/wrapper/customConverters.cpp 2014-05-23 13:05:19 +0000
+++ py/wrapper/customConverters.cpp 2014-06-17 10:18:00 +0000
@@ -149,6 +149,38 @@
}
};
+
+
+#ifdef YADE_MASK_ARBITRARY
+struct custom_mask_to_long{
+ static PyObject* convert(const mask_t& mask){
+ return PyLong_FromString(const_cast<char*>(mask.to_string().c_str()),NULL,2);
+ }
+};
+struct custom_mask_from_long{
+ custom_mask_from_long(){
+ boost::python::converter::registry::push_back(&convertible,&construct,boost::python::type_id<mask_t>());
+ }
+ static void* convertible(PyObject* obj_ptr){
+ return (PyLong_Check(obj_ptr) || PyInt_Check(obj_ptr))? obj_ptr : 0;
+ }
+ static void construct(PyObject* obj_ptr, boost::python::converter::rvalue_from_python_stage1_data* data){
+ void* storage=((boost::python::converter::rvalue_from_python_storage<mask_t>*)(data))->storage.bytes;
+ new (storage) mask_t; mask_t* mask=(mask_t*)storage;
+ if (PyInt_Check(obj_ptr)) obj_ptr = PyLong_FromLong(PyInt_AsLong(obj_ptr));
+ obj_ptr = _PyLong_Format(obj_ptr,2,0,0);
+ std::string s(PyString_AsString(obj_ptr));
+ //
+ if (s.substr(0,2).compare("0b")==0) s = s.substr(2);
+ if (s[s.length()-1]=='L') s = s.substr(0,s.length()-1);
+ // TODO?
+ *mask = mask_t(s);
+ data->convertible=storage;
+ }
+};
+#endif
+
+
BOOST_PYTHON_MODULE(_customConverters){
custom_Se3r_from_seq(); boost::python::to_python_converter<Se3r,custom_se3_to_tuple>();
@@ -167,6 +199,11 @@
//boost::python::to_python_converter<std::list<shared_ptr<Functor> >, custom_list_to_list<shared_ptr<Functor> > >();
//boost::python::to_python_converter<std::list<shared_ptr<Functor> >, custom_list_to_list<shared_ptr<Functor> > >();
+#ifdef YADE_MASK_ARBITRARY
+ custom_mask_from_long();
+ boost::python::to_python_converter<mask_t,custom_mask_to_long>();
+#endif
+
// register 2-way conversion between c++ vector and python homogeneous sequence (list/tuple) of corresponding type
#define VECTOR_SEQ_CONV(Type) custom_vector_from_seq<Type>(); boost::python::to_python_converter<std::vector<Type>, custom_vector_to_list<Type> >();
VECTOR_SEQ_CONV(int);
=== modified file 'py/ymport.py'
--- py/ymport.py 2014-05-26 07:53:31 +0000
+++ py/ymport.py 2014-06-12 16:54:31 +0000
@@ -48,7 +48,7 @@
raise RuntimeError("Please, specify a correct format output!");
return ret
-def textClumps(fileName,shift=Vector3.Zero,discretization=0,orientation=Quaternion.Identity,scale=1.0,**kw):
+def textClumps(fileName,shift=Vector3.Zero,discretization=0,orientation=Quaternion((0,1,0),0.0),scale=1.0,**kw):
"""Load clumps-members from file, insert them to the simulation.
:param str filename: file name
@@ -140,7 +140,7 @@
surf.translate(shift)
yade.pack.gtsSurface2Facets(surf,**kw)
-def gmsh(meshfile="file.mesh",shift=Vector3.Zero,scale=1.0,orientation=Quaternion.Identity,**kw):
+def gmsh(meshfile="file.mesh",shift=Vector3.Zero,scale=1.0,orientation=Quaternion((0,1,0),0.0),**kw):
""" Imports geometry from mesh file and creates facets.
:Parameters:
@@ -211,7 +211,7 @@
ret.append(utils.facet((nodelistVector3[i[1]],nodelistVector3[i[2]],nodelistVector3[i[3]]),**kw))
return ret
-def gengeoFile(fileName="file.geo",shift=Vector3.Zero,scale=1.0,orientation=Quaternion.Identity,**kw):
+def gengeoFile(fileName="file.geo",shift=Vector3.Zero,scale=1.0,orientation=Quaternion((0,1,0),0.0),**kw):
""" Imports geometry from LSMGenGeo .geo file and creates spheres.
Since 2012 the package is available in Debian/Ubuntu and known as python-demgengeo
http://packages.qa.debian.org/p/python-demgengeo.html
=== added file 'scripts/update_names.sh'
--- scripts/update_names.sh 1970-01-01 00:00:00 +0000
+++ scripts/update_names.sh 2014-06-18 17:59:20 +0000
@@ -0,0 +1,67 @@
+#!/bin/sh
+# The script modifies authors names not to get duplicated names
+# Modifies the git directory!
+
+git filter-branch --commit-filter '
+ if [ "$GIT_AUTHOR_EMAIL" = "gladk@xxxxxxxxxx" ];
+ then
+ GIT_AUTHOR_EMAIL="gladky.anton@xxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_EMAIL" = "bruno.chareyre@xxxxxxxxxxxxxxx" ]
+ then
+ GIT_AUTHOR_EMAIL="bruno.chareyre@xxxxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_EMAIL" = "christian@localhost.localdomain" ]
+ then
+ GIT_AUTHOR_NAME="Christian Jakob";
+ GIT_AUTHOR_EMAIL="jakob@xxxxxxxxxxxxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_EMAIL" = "donia.marzougui@xxxxxxxxxxxxxxx" ] || [ "$GIT_AUTHOR_NAME" = "dmarzougui" ] || [ "$GIT_AUTHOR_NAME" = "Donia" ] || [ "$GIT_AUTHOR_EMAIL" = "marzougui.donia@xxxxxxxxxx" ]
+ then
+ GIT_AUTHOR_NAME="Donia Marzougui";
+ GIT_AUTHOR_EMAIL="donia.marzougui@xxxxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Francois" ] ||[ "$GIT_AUTHOR_NAME" = "françois" ] || [ "$GIT_AUTHOR_NAME" = "Francois Kneib" ] || [ "$GIT_AUTHOR_NAME" = "François Kneib" ]
+ then
+ GIT_AUTHOR_NAME="Francois Kneib";
+ GIT_AUTHOR_EMAIL="francois.kneib@xxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Jan Stransky" ] ||[ "$GIT_AUTHOR_NAME" = "Jan Stránský" ] || [ "$GIT_AUTHOR_NAME" = "Jan Stransky" ] || [ "$GIT_AUTHOR_NAME" = "Jan Stránský" ]
+ then
+ GIT_AUTHOR_NAME="Jan Stransky";
+ GIT_AUTHOR_EMAIL="jan.stransky@xxxxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Jerome Duriez" ]
+ then
+ GIT_AUTHOR_NAME="Jerome Duriez";
+ GIT_AUTHOR_EMAIL="duriez@xxxxxxxxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Luc Scholtes" ] || [ "$GIT_AUTHOR_NAME" = "scholtes" ]
+ then
+ GIT_AUTHOR_NAME="Luc Scholtes";
+ GIT_AUTHOR_EMAIL="lscholtes63@xxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Luc Sibille" ]
+ then
+ GIT_AUTHOR_NAME="Luc Sibille";
+ GIT_AUTHOR_EMAIL="luc.sibille@xxxxxxxxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Raphael Maurin" ]
+ then
+ GIT_AUTHOR_NAME="Raphaël Maurin";
+ GIT_AUTHOR_EMAIL="raphael.maurin@xxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Raphael Maurin" ]
+ then
+ GIT_AUTHOR_NAME="Raphaël Maurin";
+ GIT_AUTHOR_EMAIL="raphael.maurin@xxxxxxxxx";
+ git commit-tree "$@";
+ elif [ "$GIT_AUTHOR_NAME" = "Kubeu" ]
+ then
+ GIT_AUTHOR_NAME="Alexander Eulitz";
+ GIT_AUTHOR_EMAIL="alexander.eulitz@xxxxxxxxxxxxxxxx";
+ git commit-tree "$@";
+ else
+ git commit-tree "$@";
+ fi' HEAD
+
