yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11382
[Branch ~yade-pkg/yade/git-trunk] Rev 3378: Merge github.com:yade/trunk into chaoUnsat
Merge authors:
Anton Gladky (gladky-anton)
Bruno Chareyre (bruno-chareyre)
Christian Jakob (jakob-ifgt)
Jan Stránský (honzik)
Klaus Thoeni (klaus.thoeni)...
------------------------------------------------------------
revno: 3378 [merge]
committer: Chao Yuan <chaoyuan2012@xxxxxxxxx>
timestamp: Thu 2013-07-04 19:06:02 +0200
message:
Merge github.com:yade/trunk into chaoUnsat
removed:
examples/clumps/adaptClumpMasses-example.py
added:
cMake/FindMetis.cmake
cMake/FindOpenBlas.cmake
lib/triangulation/FlowBoundingSphereLinSolv.hpp
lib/triangulation/FlowBoundingSphereLinSolv.ipp
lib/triangulation/PeriodicFlowLinSolv.hpp
lib/triangulation/PeriodicFlowLinSolv.ipp
scripts/checks-and-tests/checks/checkWirePM.py
modified:
CMakeLists.txt
cMake/FindCholmod.cmake
cMake/FindLoki.cmake
core/Clump.cpp
core/Clump.hpp
doc/references.bib
doc/yade-articles.bib
examples/clumps/apply-buoyancy-clumps.py
gui/CMakeLists.txt
lib/base/Math.hpp
lib/triangulation/FlowBoundingSphere.cpp
lib/triangulation/FlowBoundingSphere.hpp
lib/triangulation/FlowBoundingSphere.ipp
pkg/dem/InelastCohFrictMat.hpp
pkg/dem/InelastCohFrictPhys.cpp
pkg/dem/InelastCohFrictPhys.hpp
pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp
pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp
pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp
pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp
pkg/dem/Shop.cpp
pkg/dem/WirePM.cpp
pkg/dem/WirePM.hpp
py/CMakeLists.txt
py/wrapper/yadeWrapper.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'CMakeLists.txt'
--- CMakeLists.txt 2013-05-14 21:24:11 +0000
+++ CMakeLists.txt 2013-07-02 15:15:33 +0000
@@ -12,6 +12,8 @@
# ENABLE_OPENMP: enable OpenMP-parallelizing option (ON by default)
# ENABLE_GTS: enable GTS-option (ON by default)
# ENABLE_GL2PS: enable GL2PS-option (ON by default)
+# ENABLE_LINSOLV: enable LINSOLV-option (OFF by default)
+# ENABLE_PFVFLOW: enable PFVFLOW-option, FlowEngine (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).
@@ -89,6 +91,7 @@
OPTION(ENABLE_CGAL "Enable CGAL" ${DEFAULT})
OPTION(ENABLE_GL2PS "Enable GL2PS" ${DEFAULT})
OPTION(ENABLE_LINSOLV "Enable direct solver for the flow engines (experimental)" OFF)
+OPTION(ENABLE_PFVFLOW "Enable flow engine (experimental)" OFF)
#===========================================================
# Use Eigen3 by default
@@ -183,8 +186,6 @@
MESSAGE(STATUS "PFVFLOW depends on CGAL, attempting to turn ENABLE_CGAL ON")
SET(ENABLE_CGAL ON)
ENDIF (NOT ENABLE_CGAL)
- SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DFLOW_ENGINE")
- SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} PFVflow")
ENDIF(ENABLE_PFVFLOW)
#===========================================================
IF(ENABLE_CGAL)
@@ -201,28 +202,43 @@
SET(LINKLIBS "${LINKLIBS};${CGAL_LIBRARIES};${GMP_LIBRARIES};${GMPXX_LIBRARIES};")
MESSAGE(STATUS "Found CGAL")
SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} CGAL")
+
+ IF(ENABLE_PFVFLOW)
+ SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DFLOW_ENGINE")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} PFVflow")
+ ENDIF(ENABLE_PFVFLOW)
+
ELSE(CGAL_FOUND AND GMP_FOUND AND (NOT("${CMAKE_CXX_COMPILER} ${CMAKE_CXX_COMPILER_ARG1}" MATCHES ".*clang")))
MESSAGE(STATUS "CGAL NOT found")
SET(DISABLED_FEATS "${DISABLED_FEATS} CGAL")
SET(ENABLE_CGAL OFF)
+ IF(ENABLE_PFVFLOW)
+ SET(DISABLED_FEATS "${DISABLED_FEATS} PFVflow")
+ MESSAGE(STATUS "CGAL NOT found: PFVFLOW disabled")
+ ENDIF(ENABLE_PFVFLOW)
+
ENDIF(CGAL_FOUND AND GMP_FOUND AND (NOT("${CMAKE_CXX_COMPILER} ${CMAKE_CXX_COMPILER_ARG1}" MATCHES ".*clang")))
ELSE(ENABLE_CGAL)
SET(DISABLED_FEATS "${DISABLED_FEATS} CGAL")
ENDIF(ENABLE_CGAL)
#===========================================================
IF(ENABLE_LINSOLV)
- INCLUDE(FindCholmod)
FIND_PACKAGE(Cholmod)
- IF(CHOLMOD_FOUND)
+ FIND_PACKAGE(OpenBlas)
+ FIND_PACKAGE(Metis)
+ IF(CHOLMOD_FOUND AND OPENBLAS_FOUND AND METIS_FOUND)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CGAL_DEFINITIONS} -DLINSOLV -DFLOW_ENGINE")
- SET(LINKLIBS "${LINKLIBS};${CHOLMOD_LIBRARIES};${AMD_LIBRARY};${CAMD_LIBRARY};${COLAMD_LIBRARY};${CCOLAMD_LIBRARY};libmetis.so;libopenblas.a;")
- MESSAGE(STATUS "Found CHOLMOD")
- SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} cholmod")
- ELSE(CHOLMOD_FOUND)
+ SET(LINKLIBS "${LINKLIBS};${CHOLMOD_LIBRARIES};${AMD_LIBRARY};${CAMD_LIBRARY};${COLAMD_LIBRARY};${CCOLAMD_LIBRARY};${OPENBLAS_LIBRARY};${METIS_LIBRARY}")
+ INCLUDE_DIRECTORIES(${METIS_INCLUDE_DIR} ${CHOLMOD_INCLUDE_DIR})
+ MESSAGE(STATUS "Found Cholmod")
+ MESSAGE(STATUS "Found OpenBlas")
+ MESSAGE(STATUS "Found Metis")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} LinSolv")
+ ELSE(CHOLMOD_FOUND AND OPENBLAS_FOUND AND METIS_FOUND)
MESSAGE(STATUS "CHOLMOD NOT found, LINSOLV disabled")
- SET(DISABLED_FEATS "${DISABLED_FEATS} cholmod")
+ SET(DISABLED_FEATS "${DISABLED_FEATS} LinSolv")
SET(ENABLE_LINSOLV OFF)
- ENDIF(CHOLMOD_FOUND)
+ ENDIF(CHOLMOD_FOUND AND OPENBLAS_FOUND AND METIS_FOUND)
ENDIF(ENABLE_LINSOLV)
#===============================================
IF(ENABLE_GL2PS)
@@ -326,9 +342,7 @@
FILE(GLOB SRC_LIB_COMBINED "${CMAKE_BINARY_DIR}/lib.*.cpp")
ADD_LIBRARY(support SHARED ${SRC_LIB_COMBINED})
ELSE (CHUNKSIZE)
- ADD_LIBRARY(core SHARED ${SRC_CORE})
- ADD_LIBRARY(plugins SHARED ${SRC_PKG})
- ADD_LIBRARY(support SHARED ${SRC_LIB})
+ ADD_LIBRARY(yade SHARED ${SRC_CORE} ${SRC_PKG} ${SRC_LIB})
ENDIF (CHUNKSIZE)
#===========================================================
find_python_module(minieigen)
@@ -342,17 +356,14 @@
ADD_LIBRARY(boot SHARED ${CMAKE_CURRENT_SOURCE_DIR}/core/main/pyboot.cpp)
SET_TARGET_PROPERTIES(boot PROPERTIES PREFIX "")
-TARGET_LINK_LIBRARIES(support ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ${LINKLIBS})
-TARGET_LINK_LIBRARIES(core support ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ${LINKLIBS} -lrt)
-TARGET_LINK_LIBRARIES(boot plugins support core ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ${LINKLIBS})
+TARGET_LINK_LIBRARIES(yade ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ${LINKLIBS} -lrt)
+TARGET_LINK_LIBRARIES(boot yade)
IF(ENABLE_VTK)
- TARGET_LINK_LIBRARIES(plugins vtkHybrid)
+ TARGET_LINK_LIBRARIES(yade vtkHybrid)
ENDIF(ENABLE_VTK)
IF(ENABLE_GUI)
- TARGET_LINK_LIBRARIES(core ${GUI_LIBS})
- TARGET_LINK_LIBRARIES(plugins _GLViewer ${GUI_LIBS})
- TARGET_LINK_LIBRARIES(support ${GUI_LIBS})
+ TARGET_LINK_LIBRARIES(yade _GLViewer ${GUI_LIBS})
ENDIF(ENABLE_GUI)
#====================================
@@ -390,9 +401,7 @@
INSTALL(FILES "${CMAKE_CURRENT_SOURCE_DIR}/doc/yade-logo-note.png" DESTINATION "${YADE_DOC_PATH}/img")
INSTALL(TARGETS boot DESTINATION "${YADE_PY_PATH}/yade/")
-INSTALL(TARGETS core DESTINATION ${YADE_LIB_PATH})
-INSTALL(TARGETS plugins DESTINATION ${YADE_LIB_PATH})
-INSTALL(TARGETS support DESTINATION ${YADE_LIB_PATH})
+INSTALL(TARGETS yade DESTINATION ${YADE_LIB_PATH})
#===========================================================
MESSAGE(STATUS "===========================================================")
=== modified file 'cMake/FindCholmod.cmake'
--- cMake/FindCholmod.cmake 2012-07-03 10:54:14 +0000
+++ cMake/FindCholmod.cmake 2013-07-02 14:45:22 +0000
@@ -1,17 +1,15 @@
# - Try to find CHOLMOD
# This will define
#
-# CHOLMOD_FOUND - system has CHOLMOD
-# CHOLMOD_LIBRARIES - Link these to use CHOLMOD_FOUND
-# AMD_LIBRARY - needed by CHOLMOD
-# COLAMD_LIBRARY - needed by CHOLMOD
-# CCOLAMD_LIBRARY - needed by CHOLMOD
-# CAMD_LIBRARY - needed by CHOLMOD
+# CHOLMOD_FOUND - system has CHOLMOD
+# CHOLMOD_LIBRARIES - library to link against to use Cholmod
+# CHOLMOD_INCLUDE_DIR - where to find cholmod.h, etc.
+# AMD_LIBRARY - needed by CHOLMOD
+# COLAMD_LIBRARY - needed by CHOLMOD
+# CCOLAMD_LIBRARY - needed by CHOLMOD
+# CAMD_LIBRARY - needed by CHOLMOD
-IF (CHOLMOD_LIBRARIES)
- SET(CHOLMOD_FOUND TRUE)
-ELSE(CHOLMOD_LIBRARIES)
- FIND_LIBRARY(CHOLMOD_LIBRARIES NAMES cholmod libcholmod
+FIND_LIBRARY(CHOLMOD_LIBRARIES NAMES cholmod libcholmod
PATHS
/usr/lib
/usr/local/lib
@@ -20,15 +18,14 @@
/usr/local/lib64
/usr/lib64/CGAL
)
- IF(CHOLMOD_LIBRARIES)
- SET(CHOLMOD_FOUND TRUE)
- FIND_LIBRARY(AMD_LIBRARY NAMES amd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
- FIND_LIBRARY(CAMD_LIBRARY NAMES camd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
- FIND_LIBRARY(COLAMD_LIBRARY NAMES colamd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
- FIND_LIBRARY(CCOLAMD_LIBRARY NAMES ccolamd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
- ELSE(CHOLMOD_LIBRARIES)
- SET(CHOLMOD_FOUND FALSE)
- MESSAGE(STATUS "CHOLMOD not found.")
- ENDIF(CHOLMOD_LIBRARIES)
- MARK_AS_ADVANCED(CHOLMOD_LIBRARIES)
-ENDIF (CHOLMOD_LIBRARIES)
+
+FIND_LIBRARY(AMD_LIBRARY NAMES amd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
+FIND_LIBRARY(CAMD_LIBRARY NAMES camd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
+FIND_LIBRARY(COLAMD_LIBRARY NAMES colamd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
+FIND_LIBRARY(CCOLAMD_LIBRARY NAMES ccolamd PATHS /usr/lib /usr/local/lib /usr/lib/CGAL /usr/lib64 /usr/local/lib64 /usr/lib64/CGAL)
+
+FIND_PATH(CHOLMOD_INCLUDE_DIR cholmod.h PATH /usr/include /usr/include/suitesparse)
+
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(Cholmod DEFAULT_MSG CHOLMOD_LIBRARIES CHOLMOD_INCLUDE_DIR AMD_LIBRARY CAMD_LIBRARY COLAMD_LIBRARY CCOLAMD_LIBRARY)
+MARK_AS_ADVANCED(CHOLMOD_LIBRARIES CHOLMOD_INCLUDE_DIR AMD_LIBRARY CAMD_LIBRARY COLAMD_LIBRARY CCOLAMD_LIBRARY)
=== modified file 'cMake/FindLoki.cmake'
--- cMake/FindLoki.cmake 2013-01-29 16:16:15 +0000
+++ cMake/FindLoki.cmake 2013-06-26 18:15:55 +0000
@@ -1,5 +1,5 @@
# - Find Loki library
-# Find the native GL2PS includes and library
+#
# This module defines
# LOKI_INCLUDE_DIR, where to find loki/Typelist.h, etc.
# LOKI_LIBRARY, libraries to link against to use GL2PS.
=== added file 'cMake/FindMetis.cmake'
--- cMake/FindMetis.cmake 1970-01-01 00:00:00 +0000
+++ cMake/FindMetis.cmake 2013-06-26 18:15:55 +0000
@@ -0,0 +1,16 @@
+# - Find Metis library
+#
+# This module defines
+# METIS_INCLUDE_DIR, where to find loki/Typelist.h, etc.
+# METIS_LIBRARY, libraries to link against to use GL2PS.
+# METIS_FOUND, If false, do not try to use GL2PS.
+
+FIND_PATH(METIS_INCLUDE_DIR metis.h)
+FIND_LIBRARY(METIS_LIBRARY NAMES metis )
+
+# handle the QUIETLY and REQUIRED arguments and set LOKI_FOUND to TRUE if
+# all listed variables are TRUE
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(Metis DEFAULT_MSG METIS_INCLUDE_DIR METIS_LIBRARY)
+
+MARK_AS_ADVANCED(METIS_INCLUDE_DIR METIS_LIBRARY)
=== added file 'cMake/FindOpenBlas.cmake'
--- cMake/FindOpenBlas.cmake 1970-01-01 00:00:00 +0000
+++ cMake/FindOpenBlas.cmake 2013-06-26 18:15:55 +0000
@@ -0,0 +1,14 @@
+# - Find OpenBlas library
+#
+# This module defines
+# OPENBLAS_LIBRARY, libraries to link against to use GL2PS.
+# OPENBLAS_FOUND, If false, do not try to use GL2PS.
+
+FIND_LIBRARY(OPENBLAS_LIBRARY NAMES openblas PATHS /usr/lib/openblas-base )
+
+# handle the QUIETLY and REQUIRED arguments and set LOKI_FOUND to TRUE if
+# all listed variables are TRUE
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(OpenBlas DEFAULT_MSG OPENBLAS_LIBRARY)
+
+MARK_AS_ADVANCED(OPENBLAS_LIBRARY)
=== modified file 'core/Clump.cpp'
--- core/Clump.cpp 2013-01-31 15:44:02 +0000
+++ core/Clump.cpp 2013-07-04 14:20:19 +0000
@@ -5,7 +5,7 @@
#include<yade/core/Scene.hpp>
#include<yade/core/BodyContainer.hpp>
#include<yade/core/State.hpp>
-
+#include<yade/pkg/common/Sphere.hpp>
YADE_PLUGIN((Clump));
CREATE_LOGGER(Clump);
@@ -64,8 +64,6 @@
/*! Clump's se3 will be updated (origin at centroid and axes coincident with principal inertia axes) and subSe3 modified in such a way that members positions in world coordinates will not change.
- TODO: numerical integration of inertia based on regular space sampling with relative tolerance WRT minimum sphere radius
-
Note: velocities and angularVelocities of constituents are zeroed.
OLD DOCS (will be cleaned up):
@@ -84,7 +82,7 @@
*/
-void Clump::updateProperties(const shared_ptr<Body>& clumpBody, bool intersecting){
+void Clump::updateProperties(const shared_ptr<Body>& clumpBody){
LOG_DEBUG("Updating clump #"<<clumpBody->id<<" parameters");
const shared_ptr<State> state(clumpBody->state);
const shared_ptr<Clump> clump(YADE_PTR_CAST<Clump>(clumpBody->shape));
@@ -108,84 +106,118 @@
state->angVel=Vector3r::Zero();
return;
}
-
+ //check for intersections:
+ bool intersecting = false;
+ shared_ptr<Sphere> sph (new Sphere);
+ int Sph_Index = sph->getClassIndexStatic(); // get sphere index for checking if bodies are spheres
+ FOREACH(MemberMap::value_type& mm, clump->members){
+ const shared_ptr<Body> subBody1=Body::byId(mm.first);
+ FOREACH(MemberMap::value_type& mm, clump->members){
+ const shared_ptr<Body> subBody2=Body::byId(mm.first);
+ if ((subBody1->shape->getClassIndex() == Sph_Index) && (subBody2->shape->getClassIndex() == Sph_Index) && (subBody1!=subBody2)){//clump members should be spheres
+ Vector3r dist = subBody1->state->pos - subBody2->state->pos;
+ const Sphere* sphere1 = YADE_CAST<Sphere*> (subBody1->shape.get());
+ const Sphere* sphere2 = YADE_CAST<Sphere*> (subBody2->shape.get());
+ Real un = (sphere1->radius+sphere2->radius) - dist.norm();
+ if (un > 0.) {intersecting = true; break;}
+ }
+ }
+ }
/* quantities suffixed by
g: global (world) coordinates
s: local subBody's coordinates
c: local clump coordinates
*/
- double M=0; // mass
+ Real M=0; // mass
+ Real dens=0;//density
Vector3r Sg(0,0,0); // static moment, for getting clump's centroid
Matrix3r Ig(Matrix3r::Zero()), Ic(Matrix3r::Zero()); // tensors of inertia; is upper triangular, zeros instead of symmetric elements
-
- if(intersecting){
- LOG_WARN("Self-intersecting clumps not yet implemented, intersections will be ignored.");
- intersecting=false;
+
+ /**
+ algorithm for estimation of volumes and inertia tensor from clumps using summation/integration scheme with regular grid spacing
+ (some parts copied from woo: http://bazaar.launchpad.net/~eudoxos/woo/trunk/view/head:/pkg/dem/Clump.cpp)
+ */
+ if(intersecting){
+ //get boundaries and minimum radius of clump:
+ Real rMin=1./0.; AlignedBox3r aabb;
+ FOREACH(MemberMap::value_type& mm, clump->members){
+ const shared_ptr<Body> subBody = Body::byId(mm.first);
+ if (subBody->shape->getClassIndex() == Sph_Index){//clump member should be a sphere
+ const Sphere* sphere = YADE_CAST<Sphere*> (subBody->shape.get());
+ aabb.extend(subBody->state->pos + Vector3r::Constant(sphere->radius));
+ aabb.extend(subBody->state->pos - Vector3r::Constant(sphere->radius));
+ rMin=min(rMin,sphere->radius);
+ }
+ }
+ //get volume and inertia tensor using regular cubic cell array inside bounding box of the clump:
+ int divisor = 15; //TODO: make it choosable by users
+ Real dx = rMin/divisor; //edge length of cell
+ Real aabbMax = max(max(aabb.max().x()-aabb.min().x(),aabb.max().y()-aabb.min().y()),aabb.max().z()-aabb.min().z());
+ if (aabbMax/dx > 150) dx = aabbMax/150;//limit dx
+ Real dv = pow(dx,3); //volume of cell
+ Vector3r x; //position vector (center) of cell
+ for(x.x()=aabb.min().x()+dx/2.; x.x()<aabb.max().x(); x.x()+=dx){
+ for(x.y()=aabb.min().y()+dx/2.; x.y()<aabb.max().y(); x.y()+=dx){
+ for(x.z()=aabb.min().z()+dx/2.; x.z()<aabb.max().z(); x.z()+=dx){
+ FOREACH(MemberMap::value_type& mm, clump->members){
+ const shared_ptr<Body> subBody = Body::byId(mm.first);
+ if (subBody->shape->getClassIndex() == Sph_Index){//clump member should be a sphere
+ dens = subBody->material->density;
+ const Sphere* sphere = YADE_CAST<Sphere*> (subBody->shape.get());
+ if((x-subBody->state->pos).squaredNorm() < pow(sphere->radius,2)){
+ M += dv;
+ Sg += dv*x;
+ //inertia I = sum_i( mass_i*dist^2 + I_s) ) //steiners theorem
+ Ig += dv*( x.dot(x)*Matrix3r::Identity()-x*x.transpose()/*dist^2*/+Matrix3r(Vector3r::Constant(dv*pow(dx,2)/6.).asDiagonal())/*I_s/m = d^2: along princial axes of dv; perhaps negligible?*/);
+ }
+ }
+ }
+ }
+ }
+ }
}
-
- // begin non-intersecting loop here
if(!intersecting){
- FOREACH(MemberMap::value_type& I, clump->members){
+ FOREACH(MemberMap::value_type& mm, clump->members){
// I.first is Body::id_t, I.second is Se3r of that body
- shared_ptr<Body> subBody=Body::byId(I.first);
- State* subState=subBody->state.get();
- M+=subState->mass;
- Sg+=subState->mass*subState->pos;
- // transform from local to global coords
- Quaternionr subState_ori_conjugate=subState->ori.conjugate();
- Matrix3r Imatrix=Matrix3r::Zero(); Imatrix.diagonal()=subState->inertia;
- // TRWM3MAT(Imatrix); TRWM3QUAT(subRBP_orientation_conjugate);
- Ig+=Clump::inertiaTensorTranslate(Clump::inertiaTensorRotate(Imatrix,subState_ori_conjugate),subState->mass,-1.*subState->pos);
- //TRWM3MAT(Clump::inertiaTensorRotate(Matrix3r(subRBP->inertia),subRBP_orientation_conjugate));
+ const shared_ptr<Body> subBody=Body::byId(mm.first);
+ dens = subBody->material->density;
+ if (subBody->shape->getClassIndex() == Sph_Index){//clump member should be a sphere
+ State* subState=subBody->state.get();
+ const Sphere* sphere = YADE_CAST<Sphere*> (subBody->shape.get());
+ Real vol = (4./3.)*Mathr::PI*pow(sphere->radius,3.);
+ M+=vol;
+ Sg+=vol*subState->pos;
+ Ig+=Clump::inertiaTensorTranslate(Vector3r::Constant((2/5.)*vol*pow(sphere->radius,2)).asDiagonal(),vol,-1.*subState->pos);
+ }
}
}
- TRVAR1(M); TRWM3MAT(Ig); TRWM3VEC(Sg);
- assert(M>0);
-
- state->pos=Sg/M; // clump's centroid
+ assert(M>0); LOG_TRACE("M=\n"<<M<<"\nIg=\n"<<Ig<<"\nSg=\n"<<Sg);
+ // clump's centroid
+ state->pos=Sg/M;
// this will calculate translation only, since rotation is zero
- Matrix3r Ic_orientG=Clump::inertiaTensorTranslate(Ig, -M /* negative mass means towards centroid */, state->pos); // inertia at clump's centroid but with world orientation
- TRWM3MAT(Ic_orientG);
-
- Matrix3r R_g2c(Matrix3r::Zero()); //rotation matrix
+ Matrix3r Ic_orientG=inertiaTensorTranslate(Ig, -M /* negative mass means towards centroid */, state->pos); // inertia at clump's centroid but with world orientation
+ LOG_TRACE("Ic_orientG=\n"<<Ic_orientG);
Ic_orientG(1,0)=Ic_orientG(0,1); Ic_orientG(2,0)=Ic_orientG(0,2); Ic_orientG(2,1)=Ic_orientG(1,2); // symmetrize
- //TRWM3MAT(Ic_orientG);
- matrixEigenDecomposition(Ic_orientG,R_g2c,Ic);
- /*! @bug eigendecomposition might be wrong. see http://article.gmane.org/gmane.science.physics.yade.devel/99 for message. It is worked around below, however.
- */
- // has NaNs for identity matrix!
- TRWM3MAT(R_g2c);
-
+ Eigen::SelfAdjointEigenSolver<Matrix3r> decomposed(Ic_orientG);
+ const Matrix3r& R_g2c(decomposed.eigenvectors());
+ // has NaNs for identity matrix??
+ LOG_TRACE("R_g2c=\n"<<R_g2c);
// set quaternion from rotation matrix
state->ori=Quaternionr(R_g2c); state->ori.normalize();
- // now Ic is diagonal
- state->inertia=Ic.diagonal();
- state->mass=M;
-
-
- // this block will be removed once EigenDecomposition works for diagonal matrices
- #if 1
- if(isnan(R_g2c(0,0))||isnan(R_g2c(0,1))||isnan(R_g2c(0,2))||isnan(R_g2c(1,0))||isnan(R_g2c(1,1))||isnan(R_g2c(1,2))||isnan(R_g2c(2,0))||isnan(R_g2c(2,1))||isnan(R_g2c(2,2))){
- throw std::logic_error("Clump::updateProperties: NaNs in eigen-decomposition of inertia matrix?!");
- }
- #endif
- TRWM3VEC(state->inertia);
-
+ state->inertia=decomposed.eigenvalues();
+ state->mass=M*dens;
+
// TODO: these might be calculated from members... but complicated... - someone needs that?!
state->vel=state->angVel=Vector3r::Zero();
-
clumpBody->setAspherical(state->inertia[0]!=state->inertia[1] || state->inertia[0]!=state->inertia[2]);
// update subBodySe3s; subtract clump orientation (=apply its inverse first) to subBody's orientation
FOREACH(MemberMap::value_type& I, clump->members){
- // now, I->first is Body::id_t, I->second is Se3r of that body
shared_ptr<Body> subBody=Body::byId(I.first);
- //const shared_ptr<RigidBodyParameters>& subRBP(YADE_PTR_CAST<RigidBodyParameters>(subBody->physicalParameters));
State* subState=subBody->state.get();
I.second.orientation=state->ori.conjugate()*subState->ori;
I.second.position=state->ori.conjugate()*(subState->pos-state->pos);
}
-
}
/*! @brief Recalculates inertia tensor of a body after translation away from (default) or towards its centroid.
@@ -196,19 +228,7 @@
* @return inertia tensor in the new coordinate system; the matrix is symmetric.
*/
Matrix3r Clump::inertiaTensorTranslate(const Matrix3r& I,const Real m, const Vector3r& off){
- Real ooff=off.dot(off);
- Matrix3r I2=I;
- // TODO: replace by nicer eigen code
-
- // translation away from centroid
- /* I^c_jk=I'_jk-M*(delta_jk R.R - R_j*R_k) [http://en.wikipedia.org/wiki/Moments_of_inertia#Parallel_axes_theorem] */
- Matrix3r dI; dI<</* dIxx */ ooff-off[0]*off[0], /* dIxy */ -off[0]*off[1], /* dIxz */ -off[0]*off[2],
- /* sym */ 0., /* dIyy */ ooff-off[1]*off[1], /* dIyz */ -off[1]*off[2],
- /* sym */ 0., /* sym */ 0., /* dIzz */ ooff-off[2]*off[2];
- I2+=m*dI;
- I2(1,0)=I2(0,1); I2(2,0)=I2(0,2); I2(2,1)=I2(1,2); // symmetrize
- //TRWM3MAT(I2);
- return I2;
+ return I+m*(off.dot(off)*Matrix3r::Identity()-off*off.transpose());
}
/*! @brief Recalculate body's inertia tensor in rotated coordinates.
@@ -219,7 +239,6 @@
*/
Matrix3r Clump::inertiaTensorRotate(const Matrix3r& I,const Matrix3r& T){
/* [http://www.kwon3d.com/theory/moi/triten.html] */
- //TRWM3MAT(I); TRWM3MAT(T);
return T.transpose()*I*T;
}
@@ -233,5 +252,3 @@
Matrix3r T=rot.toRotationMatrix();
return inertiaTensorRotate(I,T);
}
-
-
=== modified file 'core/Clump.hpp'
--- core/Clump.hpp 2013-03-06 17:30:45 +0000
+++ core/Clump.hpp 2013-06-21 06:53:43 +0000
@@ -61,7 +61,7 @@
static void add(const shared_ptr<Body>& clump, const shared_ptr<Body>& subBody);
static void del(const shared_ptr<Body>& clump, const shared_ptr<Body>& subBody);
//! Recalculate physical properties of Clump.
- static void updateProperties(const shared_ptr<Body>& clump, bool intersecting);
+ static void updateProperties(const shared_ptr<Body>& clump);
//! Calculate positions and orientations of members based on relative Se3; newton pointer (if non-NULL) calls NewtonIntegrator::saveMaximaVelocity
// done as template to avoid cross-dependency between clump and newton (not necessary if all plugins are linked together)
template<class IntegratorT>
=== modified file 'doc/references.bib'
--- doc/references.bib 2013-06-12 14:35:49 +0000
+++ doc/references.bib 2013-06-25 04:00:41 +0000
@@ -587,3 +587,15 @@
keywords = "Atomic force microscopy; Equation; Plane surface; Rupture; Thermodynamics; Methodology; Gas pressure; Calculation; Solid; Interaction energy; Geometry; Plate; Separation; Prediction; Powder; Vapor; Capillary condensation; Theory; Liquid; Force",
url = "http://www.refdoc.fr/Detailnotice?idarticle=7435486";
}
+
+@article{Thoeni2013,
+ author = "K. Thoeni and C. Lambert and A. Giacomini and S.W. Sloan",
+ doi = "10.1016/j.compgeo.2012.10.014",
+ issn = "0266-352X",
+ journal = "Computers and Geotechnics",
+ pages = "158--169",
+ title = "{Discrete modelling of hexagonal wire meshes with a stochastically distorted contact model}",
+ url = "http://www.sciencedirect.com/science/article/pii/S0266352X12002121";,
+ volume = "49",
+ year = "2013"
+}
=== modified file 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2013-06-12 14:30:56 +0000
+++ doc/yade-articles.bib 2013-06-25 04:00:41 +0000
@@ -499,7 +499,71 @@
@article{Bourrier2013,
title = {Discrete modeling of granular soils reinforcement by plant roots},
journal = {Ecological Engineering},
+ author={Bourrier, F and Kneib, F and Chareyre, B and Fourcaud, T},
year = {2013},
doi = {10.1016/j.ecoleng.2013.05.002},
url = {http://dx.doi.org/10.1016/j.ecoleng.2013.05.002}
}
+
+@article{Grujicic2013,
+ title={Discrete Element Modeling and Analysis of Structural Collapse/Survivability of a Building Subjected to Improvised Explosive Device (IED) Attack},
+ author={Grujicic, M and Snipes, JS and Ramaswami, S and Yavari, R},
+ journal={Advances in Materials Science and Applications},
+ year={2013},
+ volume={2},
+ number={1},
+ pages={9--24}
+}
+
+@article{Scholtès2013,
+ title = {A \{DEM\} model for soft and hard rocks: Role of grain interlocking on strength},
+ journal = {Journal of the Mechanics and Physics of Solids},
+ volume = {61},
+ number = {2},
+ pages = {352--369},
+ year = {2013},
+ doi = {10.1016/j.jmps.2012.10.005},
+ url = {http://www.sciencedirect.com/science/article/pii/S0022509612002268},
+ author = {Luc Scholtès and Frédéric-Victor Donzé}
+}
+
+@article{Boon2012,
+ title = {A new contact detection algorithm for three-dimensional non-spherical particles},
+ journal = {Powder Technology},
+ year = {2013},
+ doi = {10.1016/j.powtec.2012.12.040},
+ url = {http://www.sciencedirect.com/science/article/pii/S003259101200839X},
+ author = {C.W. Boon and G.T. Houlsby and S. Utili}
+}
+
+@article{Favier2012,
+ title={Rigid obstacle impacted by a supercritical cohesive granular flow using a 3D discrete element model},
+ author={Favier, L. and Daudon, D. and Donzé, F.V.},
+ journal={Cold Regions Science and Technology},
+ year={2012},
+ volume = {85},
+ pages = {232--241},
+ url = {http://dx.doi.org/10.1016/j.coldregions.2012.09.010}
+}
+
+@Article{Catalano2013a,
+ title = {Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects},
+ author = {E. Catalano and B. Chareyre and E. Barthélémy},
+ journal = {International Journal for Numerical and Analytical Methods in Geomechanics},
+ year = {2013},
+ doi = {10.1002/nag.2198},
+ note = {in press, arxiv version available},
+ url = {http://arxiv.org/pdf/1304.4895.pdf}
+}
+
+@article{Thoeni2013,
+ author = "K. Thoeni and C. Lambert and A. Giacomini and S.W. Sloan",
+ doi = "10.1016/j.compgeo.2012.10.014",
+ issn = "0266-352X",
+ journal = "Computers and Geotechnics",
+ pages = "158--169",
+ title = "{Discrete modelling of hexagonal wire meshes with a stochastically distorted contact model}",
+ url = "http://www.sciencedirect.com/science/article/pii/S0266352X12002121";,
+ volume = "49",
+ year = "2013"
+}
\ No newline at end of file
=== removed file 'examples/clumps/adaptClumpMasses-example.py'
--- examples/clumps/adaptClumpMasses-example.py 2013-04-17 07:39:24 +0000
+++ examples/clumps/adaptClumpMasses-example.py 1970-01-01 00:00:00 +0000
@@ -1,54 +0,0 @@
-#!/usr/bin/python
-# -*- coding: utf-8 -*-
-
-'''This example shows usage of adaptClumpMasses().'''
-
-#define material for all bodies:
-id_Mat=O.materials.append(FrictMat(young=1e7,poisson=0.3,density=1000,frictionAngle=1))
-Mat=O.materials[id_Mat]
-
-#define engines:
-O.engines=[
- ForceResetter(),
- InsertionSortCollider([Bo1_Sphere_Aabb()]),
- InteractionLoop(
- [Ig2_Sphere_Sphere_ScGeom()],
- [Ip2_FrictMat_FrictMat_FrictPhys()],
- [Law2_ScGeom_FrictPhys_CundallStrack()]
- ),
- NewtonIntegrator(damping=0.7,gravity=[0,0,-10])
-]
-
-#append spheres:
-sphereList = []
-sphereList.append(O.bodies.append(sphere((0,0,0),radius=.5,material=Mat)))
-sphereList.append(O.bodies.append(sphere((.5,0,0),radius=.25,material=Mat)))
-sphereList.append(O.bodies.append(sphere((0,.3,0),radius=.5,material=Mat)))
-sphereList.append(O.bodies.append(sphere((0,0,.3),radius=.5,material=Mat)))
-
-#clump spheres:
-idClump=O.bodies.clump(sphereList)
-
-from yade import qt
-qt.Controller()
-qt.View()
-
-def printClumpMasses():
- for b in O.bodies:
- if b.isClump:
- print 'Clump ',b.id,' with mass ',b.state.mass,' has following members:'
- keys = b.shape.members.keys()
- for ii in range(0,len(keys)):
- print '- Body ',keys[ii],' with mass ',O.bodies[keys[ii]].state.mass
-
-print 'Mass info BEFORE adaptClumpMasses() is called: ----------------'
-printClumpMasses()
-
-massInfo = O.bodies.adaptClumpMasses([],10000000) #give an empty list [] if no body should be excluded, number of grid points is set to 10000000
-
-#NOTE: the higher the number of grid points, the more precisely is the result,
-# but high numbers of clumps + high number of grid points = a looooong coffee break ...
-
-print 'Mass info AFTER adaptClumpMasses() is called: -----------------'
-printClumpMasses()
-print 'output of adaptClumpMasses(): ',massInfo
=== modified file 'examples/clumps/apply-buoyancy-clumps.py'
--- examples/clumps/apply-buoyancy-clumps.py 2013-06-15 08:19:25 +0000
+++ examples/clumps/apply-buoyancy-clumps.py 2013-06-24 13:23:34 +0000
@@ -78,7 +78,7 @@
pos = O.bodies[keys[ii]].state.pos
h_low = min(h_low,pos[2]-O.bodies[keys[ii]].shape.radius)
h_high = max(h_high,pos[2]+O.bodies[keys[ii]].shape.radius)
- rad = ( 3*startMass[b.id]/(4*pi*O.bodies[keys[ii]].mat.density) )**(1./3.) #get radius from startMass
+ rad = ( 3*startMass[b.id]/(4*pi*O.bodies[keys[0]].mat.density) )**(1./3.) #get radius from startMass
else:
continue
if water_height > h_low:
=== modified file 'gui/CMakeLists.txt'
--- gui/CMakeLists.txt 2012-11-01 21:42:22 +0000
+++ gui/CMakeLists.txt 2013-07-01 18:57:33 +0000
@@ -9,7 +9,10 @@
ADD_LIBRARY(_GLViewer SHARED ${_GLViewer_SOURCE_FILES} ${_GLViewer_MOC_OUTFILES})
SET_TARGET_PROPERTIES(_GLViewer PROPERTIES PREFIX "")
- TARGET_LINK_LIBRARIES(_GLViewer support core ${GLUT_LIBRARY} ${OPENGL_LIBRARY} ${QT_LIBRARIES} ${QGLVIEWER_LIBRARIES} ${Boost_LIBRARIES} ${PYTHON_LIBRARIES})
+ TARGET_LINK_LIBRARIES(_GLViewer ${GLUT_LIBRARY} ${OPENGL_LIBRARY} ${QT_LIBRARIES} ${QGLVIEWER_LIBRARIES} ${Boost_LIBRARIES} ${PYTHON_LIBRARIES})
+ IF(GL2PS_FOUND AND ENABLE_GL2PS)
+ TARGET_LINK_LIBRARIES(_GLViewer ${GL2PS_LIBRARIES})
+ ENDIF(GL2PS_FOUND AND ENABLE_GL2PS)
INSTALL(TARGETS _GLViewer DESTINATION ${YADE_PY_PATH}/yade/qt)
FILE(GLOB filesPYQT "${CMAKE_CURRENT_SOURCE_DIR}/qt4/*.py")
=== modified file 'lib/base/Math.hpp'
--- lib/base/Math.hpp 2013-03-07 18:28:01 +0000
+++ lib/base/Math.hpp 2013-06-20 14:56:35 +0000
@@ -55,6 +55,8 @@
typedef Eigen::Quaternion<Real> Quaternionr;
typedef Eigen::AngleAxis<Real> AngleAxisr;
+typedef Eigen::AlignedBox<Real,2> AlignedBox2r;
+typedef Eigen::AlignedBox<Real,3> AlignedBox3r;
using Eigen::AngleAxis; using Eigen::Quaternion;
// in some cases, we want to initialize types that have no default constructor (OpenMPAccumulator, for instance)
=== modified file 'lib/triangulation/FlowBoundingSphere.cpp'
--- lib/triangulation/FlowBoundingSphere.cpp 2012-02-14 16:58:21 +0000
+++ lib/triangulation/FlowBoundingSphere.cpp 2013-06-26 18:15:55 +0000
@@ -6,7 +6,7 @@
*************************************************************************/
#ifdef FLOW_ENGINE
-#include "FlowBoundingSphere.hpp"
+#include "yade/lib/triangulation/FlowBoundingSphere.hpp"
namespace CGT {
Vecteur PeriodicCellInfo::gradP;
=== modified file 'lib/triangulation/FlowBoundingSphere.hpp'
--- lib/triangulation/FlowBoundingSphere.hpp 2013-05-03 16:45:52 +0000
+++ lib/triangulation/FlowBoundingSphere.hpp 2013-06-26 18:15:55 +0000
@@ -166,11 +166,11 @@
} //namespace CGT
#ifdef LINSOLV
-#include "FlowBoundingSphereLinSolv.hpp"
+#include "yade/lib/triangulation/FlowBoundingSphereLinSolv.hpp"
#endif
/// _____ Template Implementation ____
-#include "FlowBoundingSphere.ipp"
+#include "yade/lib/triangulation/FlowBoundingSphere.ipp"
#endif //FLOW_ENGINE
=== modified file 'lib/triangulation/FlowBoundingSphere.ipp'
--- lib/triangulation/FlowBoundingSphere.ipp 2013-06-07 07:47:59 +0000
+++ lib/triangulation/FlowBoundingSphere.ipp 2013-07-02 14:45:47 +0000
@@ -787,8 +787,8 @@
if(permeability_map){
Cell_handle c = cell;
- cell->info().s = cell->info().s + k*distance/fluidArea*Volume_Pore_VoronoiFraction (c,j);
- volume_sub_pore += Volume_Pore_VoronoiFraction (c,j);}
+ cell->info().s = cell->info().s + k*distance/fluidArea*this->Volume_Pore_VoronoiFraction (c,j);
+ volume_sub_pore += this->Volume_Pore_VoronoiFraction (c,j);}
}
}
@@ -982,8 +982,8 @@
RTriangulation& Tri = T[currentTes].Triangulation();
if (Tri.is_infinite(cell->neighbor(j))) return 0;
- double Vpore = Volume_Pore_VoronoiFraction(cell, j);
- double Ssolid = Surface_Solid_Pore(cell, j, SLIP_ON_LATERALS);
+ double Vpore = this->Volume_Pore_VoronoiFraction(cell, j);
+ double Ssolid = this->Surface_Solid_Pore(cell, j, SLIP_ON_LATERALS);
//handle symmetry (tested ok)
if (SLIP_ON_LATERALS && fictious_vertex>0) {
=== added file 'lib/triangulation/FlowBoundingSphereLinSolv.hpp'
--- lib/triangulation/FlowBoundingSphereLinSolv.hpp 1970-01-01 00:00:00 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.hpp 2013-06-27 11:29:17 +0000
@@ -0,0 +1,192 @@
+/*************************************************************************
+* Copyright (C) 2010 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 FLOW_ENGINE
+#pragma once
+
+#define EIGENSPARSE_LIB //comment this if CHOLMOD is not available
+// #define TAUCS_LIB //comment this if TAUCS lib is not available, it will disable PARDISO lib as well
+
+#ifdef EIGENSPARSE_LIB
+ #include <eigen3/Eigen/Sparse>
+ #include <eigen3/Eigen/SparseCore>
+ #include <eigen3/Eigen/CholmodSupport>
+#endif
+#ifdef TAUCS_LIB
+#define TAUCS_CORE_DOUBLE
+#include <complex> //THIS ONE MUST ABSOLUTELY BE INCLUDED BEFORE TAUCS.H!
+#include <stdlib.h>
+#include <float.h>
+extern "C" {
+#include "taucs.h"
+}
+#endif
+
+#include "FlowBoundingSphere.hpp"
+
+///_____ Declaration ____
+
+using namespace std;
+
+namespace CGT {
+
+template<class FlowType>
+class FlowBoundingSphereLinSolv : public FlowType
+{
+public:
+ DECLARE_TESSELATION_TYPES(FlowType)
+ typedef typename FlowType::Tesselation Tesselation;
+ using FlowType::useSolver;
+ using FlowType::T;
+ using FlowType::currentTes;
+ using FlowType::boundary;
+ using FlowType::y_min_id;
+ using FlowType::y_max_id;
+ using FlowType::DEBUG_OUT;
+ using FlowType::TOLERANCE;
+ using FlowType::RELAX;
+ using FlowType::fluidBulkModulus;
+ using FlowType::reApplyBoundaryConditions;
+ using FlowType::pressureChanged;
+ using FlowType::computedOnce;
+
+ //! TAUCS DECs
+ vector<Finite_cells_iterator> orderedCells;
+ bool isLinearSystemSet;
+ bool isFullLinearSystemGSSet;
+ bool areCellsOrdered;//true when orderedCells is filled, turn it false after retriangulation
+
+ #ifdef EIGENSPARSE_LIB
+ //Eigen's sparse matrix and solver
+ Eigen::SparseMatrix<double> A;
+ typedef Eigen::Triplet<double> ETriplet;
+ std::vector<ETriplet> tripletList;//The list of non-zero components in Eigen sparse matrix
+ Eigen::CholmodDecomposition<Eigen::SparseMatrix<double>, Eigen::Lower > eSolver;
+ bool factorizedEigenSolver;
+ void exportMatrix(const char* filename) {ofstream f; f.open(filename); f<<A; f.close();};
+ void exportTriplets(const char* filename) {ofstream f; f.open(filename);
+ for (int k=0; k<A.outerSize(); ++k)
+ for (Eigen::SparseMatrix<double>::InnerIterator it(A,k); it; ++it) f<< it.row()<<" "<< it.col()<<" "<<it.value()<<endl; f.close();};
+ //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 both thread numbers independently
+ int numFactorizeThreads;
+ int numSolveThreads;
+ #endif
+
+ #ifdef TAUCS_LIB
+ taucs_ccs_matrix SystemMatrix;
+ taucs_ccs_matrix* T_A;
+ taucs_ccs_matrix* Fccs;
+ void* F;//The taucs factor
+ #endif
+
+
+ int T_nnz;
+ int ncols;
+ int T_size;
+
+ double pTime1, pTime2;
+ int pTimeInt, pTime1N, pTime2N;
+
+ double ZERO;
+ vector<double> T_an;//(size*5);
+ vector<int> T_jn;//(size+1);
+ vector<int> T_ia;//(size*5);
+ vector<double> T_f;//(size); // right-hand size vector object
+ vector<Cell_handle> T_cells;//(size)
+ int T_index;
+
+ vector<double> T_b;
+ vector<double> T_bv;
+ vector <double> T_x, P_x;
+ vector <double> bodv;
+ vector <double> xodv;
+ int* perm;
+ int* invperm;
+ bool pardisoInitialized;
+ //! END TAUCS DECs
+
+
+ //! Pardiso
+ int* ia;
+ int* ja;
+ double* a;
+ int nnz;
+ int mtype; /* Real symmetric positive def. matrix */
+ double* b;
+ double* x;// the unknown vector to solve Ax=b
+ int nrhs; /* Number of right hand sides. */
+ void *pt[64];
+ int iparm[64];
+ double dparm[64];
+ int maxfct, mnum, phase, error, msglvl, solver;
+ int num_procs;
+ char *var;
+ int i;
+ double ddum; /* Double dummy */
+ int idum; /* Integer dummy. */
+ //! end pardiso
+
+ /// EXTERNAL_GS part
+ vector<vector<double> > fullAvalues;//contains Kij's and 1/(sum Kij) in 5th value (for use in GuaussSeidel)
+ vector<vector<double*> > fullAcolumns;//contains columns numbers
+ vector<double> gsP;//a vector of pressures
+ vector<double> gsdV;//a vector of dV
+ vector<double> gsB;//a vector of dV
+
+public:
+ virtual ~FlowBoundingSphereLinSolv();
+ FlowBoundingSphereLinSolv();
+
+ ///Linear system solve
+ virtual int SetLinearSystem(Real dt);
+ void VectorizedGaussSeidel(Real dt);
+ virtual int SetLinearSystemFullGS(Real dt);
+
+ int TaucsSolveTest();
+ int TaucsSolve(Real dt);
+ int PardisoSolveTest();
+ int PardisoSolve(Real dt);
+ int eigenSolve(Real dt);
+
+ void CopyGsToCells();
+ void CopyCellsToGs(Real dt);
+
+ void CopyLinToCells();
+ void CopyCellsToLin(Real dt);
+ void swap_fwd (double* v, int i);
+ void swap_fwd (int* v, int i);
+ void sort_v(int k1, int k2, int* is, double* ds);
+
+ virtual void GaussSeidel (Real dt) {
+ switch (useSolver) {
+ case 0:
+ VectorizedGaussSeidel(dt);
+ break;
+ case 1:
+ TaucsSolve(dt);
+ break;
+ case 2:
+ PardisoSolve(dt);
+ break;
+ case 3:
+ eigenSolve(dt);
+ break;
+ }
+ computedOnce=true;
+ }
+ virtual void ResetNetwork();
+};
+
+} //namespace CGT
+
+
+///_____ Implementation ____
+
+#include "FlowBoundingSphereLinSolv.ipp"
+
+#endif
=== added file 'lib/triangulation/FlowBoundingSphereLinSolv.ipp'
--- lib/triangulation/FlowBoundingSphereLinSolv.ipp 1970-01-01 00:00:00 +0000
+++ lib/triangulation/FlowBoundingSphereLinSolv.ipp 2013-06-27 11:29:17 +0000
@@ -0,0 +1,995 @@
+/*************************************************************************
+* Copyright (C) 2010 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 FLOW_ENGINE
+
+// #define XVIEW
+#include "FlowBoundingSphereLinSolv.hpp"//include after #define XVIEW
+// #include "def_types.h"
+// #include "def_flow_types.h"
+#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#include <CGAL/Width_3.h>
+#include <iostream>
+#include <fstream>
+#include <new>
+#include <utility>
+#include "vector"
+#include <assert.h>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <omp.h>
+
+
+#ifdef XVIEW
+//#include "Vue3D.h" //FIXME implicit dependencies will look for this class (out of tree) even ifndef XVIEW
+#endif
+
+#ifdef YADE_OPENMP
+// #define GS_OPEN_MP //It should never be defined if Yade is not using openmp
+#endif
+
+// #define PARDISO //comment this if pardiso lib is not available
+
+#ifdef EIGENSPARSE_LIB
+extern "C" { void openblas_set_num_threads(int num_threads); }
+#endif
+
+namespace CGT
+{
+
+using namespace std;
+
+#ifdef PARDISO
+#ifdef AIX
+#define F77_FUNC(func) func
+#else
+#define F77_FUNC(func) func ## _
+#endif
+/* PARDISO prototype. */
+extern "C" int F77_FUNC(pardisoinit)
+ (void *, int *, int *, int *, double *, int *);
+
+extern "C" int F77_FUNC(pardiso)
+ (void *, int *, int *, int *, int *, int *,
+ double *, int *, int *, int *, int *, int *,
+ int *, double *, double *, int *, double *);
+#endif
+
+#ifdef XVIEW
+Vue3D Vue1;
+#endif
+template<class FlowType>
+FlowBoundingSphereLinSolv<FlowType>::~FlowBoundingSphereLinSolv()
+{
+ #ifdef TAUCS_LIB
+ if (Fccs) taucs_ccs_free(Fccs);
+ #endif
+}
+template<class FlowType>
+FlowBoundingSphereLinSolv<FlowType>::FlowBoundingSphereLinSolv(): FlowType() {
+ useSolver=0;
+ isLinearSystemSet=0;
+ isFullLinearSystemGSSet=0;
+ areCellsOrdered=0;//true when orderedCells is filled, turn it false after retriangulation
+ ZERO=0;
+ #ifdef TAUCS_LIB
+ T_A = &SystemMatrix;
+ F = NULL;//The taucs factor
+ Fccs = NULL;//The taucs factor in CCS format
+ #endif
+ pardisoInitialized=false;
+ pTimeInt=0;pTime1N=0;pTime2N=0;
+ pTime1=0;pTime2=0;
+ #ifdef EIGENSPARSE_LIB
+ factorizedEigenSolver=false;
+ numFactorizeThreads=1;
+ numSolveThreads=1;
+ #endif
+}
+
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::swap_fwd (double* v, int i) {double temp = v[i]; v[i]=v[i+1]; v[i+1]=temp;}
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::swap_fwd (int* v, int i) {double temp = v[i]; v[i]=v[i+1]; v[i+1]=temp;}
+
+//spatial sort traits to use with a pair of CGAL::sphere pointers and integer.
+//template<class _Triangulation>
+// usage : spatial_sort(pointsPtrs.begin(),pointsPtrs.end(), RTraits_for_spatial_sort()/*, CGT::RTriangulation::Weighted_point*/);
+#include <utility>
+template<class Triangulation>
+struct CellTraits_for_spatial_sort : public Triangulation::Geom_traits {
+ typedef typename Triangulation::Geom_traits Gt;
+ typedef const typename Triangulation::Finite_cells_iterator Point_3;
+ struct Less_x_3 {bool operator()(const Point_3& p,const Point_3& q) const {
+ return typename Gt::Less_x_3()(p->info(),q->info());}
+ };
+ struct Less_y_3 {bool operator()(const Point_3& p,const Point_3& q) const {
+ return typename Gt::Less_y_3()(p->info(),q->info());}
+ };
+ struct Less_z_3 {bool operator()(const Point_3& p,const Point_3& q) const {
+ return typename Gt::Less_z_3()(p->info(),q->info());}
+ };
+ Less_x_3 less_x_3_object() const {return Less_x_3();}
+ Less_y_3 less_y_3_object() const {return Less_y_3();}
+ Less_z_3 less_z_3_object() const {return Less_z_3();}
+};
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::ResetNetwork() {
+ FlowType::ResetNetwork();
+ isLinearSystemSet=false;
+ isFullLinearSystemGSSet=false;
+ areCellsOrdered=false;
+#ifdef TAUCS_LIB
+ if (F) taucs_supernodal_factor_free(F); F=NULL;
+ if (Fccs) taucs_ccs_free(Fccs); Fccs=NULL;
+#endif
+#ifdef EIGENSPARSE_LIB
+ factorizedEigenSolver=false;
+#endif
+#ifdef PARDISO
+ if (pardisoInitialized) {
+ phase = -1;
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &ncols, &ddum, NULL, NULL, &idum, &nrhs,
+ iparm, &msglvl, &ddum, &ddum, &error, dparm);
+ pardisoInitialized=false;
+ }
+#endif
+}
+template<class FlowType>
+int FlowBoundingSphereLinSolv<FlowType>::SetLinearSystem(Real dt)
+{
+
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ int n_cells=Tri.number_of_finite_cells();
+ vector<int> clen;
+ vector<int> is;
+ vector<int> js;
+ vector<double> vs;
+ if (!areCellsOrdered) {
+ T_nnz=0;
+ ncols=0;
+ ///Ordered cells
+ orderedCells.clear();
+ const Finite_cells_iterator cell_end = Tri.finite_cells_end();
+ for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ orderedCells.push_back(cell);
+ if (!cell->info().Pcondition) ++ncols;}
+ spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
+
+ T_cells.clear();
+ T_index=0;
+ isLinearSystemSet=false;
+ areCellsOrdered=true;
+ }
+ if (!isLinearSystemSet) {
+ #ifdef TAUCS_LIB
+ if (Fccs) taucs_ccs_free(Fccs);//delete the old factor
+ #endif
+ int n = 3*(ncols+1);//number of non-zero in triangular matrix
+ is.resize(n);
+ js.resize(n);
+ vs.resize(n);
+ T_x.resize(ncols);
+ T_b.resize(ncols);
+ T_bv.resize(ncols);
+ bodv.resize(ncols);
+ xodv.resize(ncols);
+// gsB.resize(ncols+1);
+ T_cells.resize(ncols+1);
+ T_nnz=0;}
+ for (int kk=0; kk<ncols;kk++) T_b[kk]=0;
+ ///Ordered cells
+ int index=0, nIndex=0; Cell_handle neighbour_cell;
+ for (int i=0; i<n_cells; i++)
+ {
+ Finite_cells_iterator& cell = orderedCells[i];
+ ///Non-ordered cells
+// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ if (!cell->info().Pcondition) {
+ index=cell->info().index;
+ if (index==0) {
+ T_cells[++T_index]=cell;
+ cell->info().index=index=T_index;
+ }
+ if (!isLinearSystemSet) {
+ //Add diagonal term
+ is[T_nnz] = index;
+ js[T_nnz] = index;
+ vs[T_nnz] = (cell->info().k_norm())[0]+ (cell->info().k_norm())[1]+ (cell->info().k_norm())[2]+ (cell->info().k_norm())[3];
+ if (fluidBulkModulus>0) vs[T_nnz] += (1.f/(dt*fluidBulkModulus*cell->info().invVoidVolume()));
+ ++T_nnz;
+ }
+ for (int j=0; j<4; j++) {
+ neighbour_cell = cell->neighbor(j);
+ nIndex=neighbour_cell->info().index;
+ if (Tri.is_infinite(neighbour_cell)) continue;
+ if (!isLinearSystemSet && !neighbour_cell->info().Pcondition) {
+ if (nIndex==0) {
+ T_cells[++T_index]=neighbour_cell;
+ neighbour_cell->info().index=nIndex=T_index;
+ } else if (index > nIndex) {
+ is[T_nnz] = index;
+ js[T_nnz] = nIndex;
+ vs[T_nnz] = - (cell->info().k_norm())[j];
+ T_nnz++;
+ }
+ } else if (neighbour_cell->info().Pcondition) {
+ //ADD TO b, FIXME : use updated volume change
+ T_b[index-1]+=cell->info().k_norm()[j]*neighbour_cell->info().p();
+ }
+ }
+ }
+ }
+ if (!isLinearSystemSet) {
+ if (useSolver==1 || useSolver==2){
+ #ifdef TAUCS_LIB
+ clen.resize(ncols+1);
+ T_jn.resize(ncols+1);
+ T_A->colptr = &T_jn[0];
+ T_ia.resize(T_nnz);
+ T_A->rowind = &T_ia[0];
+ T_A->flags = (TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
+ T_an.resize(T_nnz);
+ T_A->values.d = &T_an[0];
+ T_A->n = ncols;
+ T_A->m = ncols;
+ int i,j,k;
+ for (j=0; j<ncols; j++) clen[j] = 0;
+ for (k=0; k<T_nnz; k++) {
+ i = is[k]-1; /* make it 1-based */
+ j = js[k]-1; /* make it 1-based */
+ (clen[j])++;
+ }
+ /* now compute column pointers */
+ k = 0;
+ for (j=0; j<ncols; j++) {
+ int tmp;
+ tmp = clen[j];
+ clen[j] = (T_A->colptr[j]) = k;
+ k += tmp;
+ }
+ clen[ncols] = (T_A->colptr[ncols]) = k;
+
+ /* now read matrix into data structure */
+ for (k=0; k<T_nnz; k++) {
+ i = is[k] - 1; /* make it 1-based */
+ j = js[k] - 1; /* make it 1-based */
+ assert(i < ncols);
+ assert(j < ncols);
+ (T_A->taucs_values)[clen[j]] = vs[k];
+ (T_A->rowind)[clen[j]] = i;
+ clen[j] ++;
+ // cerr<<"i="<< i <<" j="<< j<<" v="<<vs[k]<<" clen[j]="<<clen[j]-1<<endl;
+ }
+ #endif //TAUCS_LIB
+ #ifdef EIGENSPARSE_LIB
+ } else if (useSolver==3){
+ tripletList.clear(); tripletList.resize(T_nnz);
+ for(int k=0;k<T_nnz;k++) {
+// if (is[k]<js[k]) cerr<<"not the good relation"<<endl;
+// else cerr<< "comp " <<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
+ tripletList[k]=ETriplet(is[k]-1,js[k]-1,vs[k]);
+ }
+ A.resize(ncols,ncols);
+ A.setFromTriplets(tripletList.begin(), tripletList.end());
+ #endif
+ }
+ isLinearSystemSet=true;
+ }
+ return ncols;
+}
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::CopyGsToCells() {for (int ii=1; ii<=ncols; ii++) T_cells[ii]->info().p()=gsP[ii];}
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::CopyCellsToGs (Real dt)
+{
+ for (int ii=1; ii<=ncols; ii++){
+ gsP[ii]=T_cells[ii]->info().p();
+ gsdV[ii]= T_cells[ii]->info().dv();
+ if (fluidBulkModulus>0) { gsdV[ii] -= T_cells[ii]->info().p()/(fluidBulkModulus*dt*T_cells[ii]->info().invVoidVolume());}
+ }
+}
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::CopyLinToCells() {for (int ii=1; ii<=ncols; ii++) {T_cells[ii]->info().p()=T_x[ii-1];} }
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::CopyCellsToLin (Real dt)
+{
+ for (int ii=1; ii<=ncols; ii++) {
+ T_bv[ii-1]=T_b[ii-1]-T_cells[ii]->info().dv();
+ if (fluidBulkModulus>0) T_bv[ii-1] += T_cells[ii]->info().p()/(fluidBulkModulus*dt*T_cells[ii]->info().invVoidVolume());}
+}
+
+/// For Gauss Seidel, we need the full matrix
+template<class FlowType>
+// int FlowBoundingSphereLinSolv<FlowType>::SetLinearSystemFullGS()
+int FlowBoundingSphereLinSolv<FlowType>::SetLinearSystemFullGS(Real dt)
+{
+ //WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ int n_cells=Tri.number_of_finite_cells();
+ if (!areCellsOrdered) {
+ T_cells.clear();
+ T_index=0;
+ T_nnz=0;
+ ncols=0;
+ const Finite_cells_iterator cell_end = Tri.finite_cells_end();
+ orderedCells.clear();
+
+ for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ orderedCells.push_back(cell);
+ if (!cell->info().Pcondition) ++ncols;
+ }
+ //FIXME: does it really help? test by commenting this "sorting" line
+ spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort<RTriangulation>());
+
+// double P_zero=0;
+// if (y_min_id>=0 and y_max_id>y_min_id) P_zero = abs((boundary(y_min_id).value-boundary(y_max_id).value)/2);
+ gsP.resize(ncols+1);
+// _gsP.resize(ncols+1);
+ gsB.resize(ncols+1);
+ T_b.resize(ncols+1);
+ gsdV.resize(ncols+1);
+ fullAcolumns.resize(ncols+1);
+ fullAvalues.resize(ncols+1);
+ T_cells.resize(ncols+1);
+ for (int k=0; k<=ncols;k++) {
+ fullAcolumns[k].resize(4);
+ fullAvalues[k].resize(5);
+ gsdV[k]=0;
+// gsP[k]=P_zero;
+ }
+// _gsP[0]= &ZERO;
+ gsP[0]=0;
+ areCellsOrdered=true;
+ isFullLinearSystemGSSet=false;
+ }
+ for (int k=0; k<=ncols;k++) gsB[k]=0;
+
+ ///we build the full matrix + RHS here, else only the RHS in the other loop
+ if (!isFullLinearSystemGSSet)
+ ///Ordered cells
+ for (int i=0; i<n_cells; i++)
+ {
+ Finite_cells_iterator& cell = orderedCells[i];
+ ///Non-ordered cells
+// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ if (!cell->info().Pcondition) {
+ if (cell->info().index==0) {
+ T_cells[++T_index]=cell;
+ cell->info().index=T_index;
+ }
+ gsP[cell->info().index]=cell->info().pression;
+ //Add diagonal term
+ double num = (cell->info().k_norm())[0]+ (cell->info().k_norm())[1]+ (cell->info().k_norm())[2]+ (cell->info().k_norm())[3];
+ if (fluidBulkModulus>0) num += (1.f/(dt*fluidBulkModulus*cell->info().invVoidVolume()));
+ fullAvalues[cell->info().index][4] = 1.f/num;
+ ++T_nnz;
+
+ for (int j=0; j<4; j++) {
+ Cell_handle neighbour_cell = cell->neighbor(j);
+ if (Tri.is_infinite(neighbour_cell)) {
+ fullAvalues[cell->info().index][j] = 0;
+ fullAcolumns[cell->info().index][j] = &gsP[0];
+ continue;}
+ if (!neighbour_cell->info().Pcondition) {
+ if (neighbour_cell->info().index==0) {
+ T_cells[++T_index]=neighbour_cell;
+ neighbour_cell->info().index=T_index;
+ }
+ ++T_nnz;
+ fullAvalues[cell->info().index][j] = (cell->info().k_norm())[j];
+ fullAcolumns[cell->info().index][j] = &gsP[neighbour_cell->info().index];
+ } else {
+ fullAvalues[cell->info().index][j] = 0;
+ fullAcolumns[cell->info().index][j] = &gsP[0];
+ gsB[cell->info().index]+=cell->info().k_norm()[j]*neighbour_cell->info().p();
+ }
+ }
+ }
+ }
+ ///define only the new RHS, accouting for new imposed pressures
+ else
+ ///Ordered cells
+ for (int i=0; i<n_cells; i++)
+ {
+ Finite_cells_iterator& cell = orderedCells[i];
+ ///Non-ordered cells
+// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ if (!cell->info().Pcondition) for (int j=0; j<4; j++) {
+ Cell_handle neighbour_cell = cell->neighbor(j);
+ if (!Tri.is_infinite(neighbour_cell) && neighbour_cell->info().Pcondition)
+ gsB[cell->info().index]+=cell->info().k_norm()[j]*neighbour_cell->info().p();
+ }
+ }
+ isFullLinearSystemGSSet=true;
+ return ncols;
+}
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::VectorizedGaussSeidel(Real dt)
+{
+// cout<<"VectorizedGaussSeidel"<<endl;
+ if (!isFullLinearSystemGSSet || (isFullLinearSystemGSSet && reApplyBoundaryConditions())) SetLinearSystemFullGS(dt);
+ CopyCellsToGs(dt);
+
+ int j = 0;
+ double dp_max, p_max, sum_p, p_moy, dp_moy, sum_dp;
+ double tolerance = TOLERANCE;
+ double relax = RELAX;
+
+#ifdef GS_OPEN_MP
+ const int num_threads=1;
+ omp_set_num_threads(num_threads);
+ vector<Real> t_sum_p, t_dp_max, t_sum_dp, t_p_max;
+ t_sum_dp.resize(num_threads);
+ t_dp_max.resize(num_threads);
+ t_p_max.resize(num_threads);
+ t_sum_p.resize(num_threads);
+#endif
+ int j2=-1;
+ dp_max = 0; p_max = 0; p_moy=0; dp_moy=0; sum_p=0; sum_dp=0;
+ do {
+ if (++j2>=10) j2=0;//compute max/mean only each 10 iterations
+ if (j2==0) {
+ dp_max = 0; p_max = 0; p_moy=0; dp_moy=0; sum_p=0; sum_dp=0;
+#ifdef GS_OPEN_MP
+ for (int ii=0;ii<num_threads;ii++) t_p_max[ii] =0;
+ for (int ii=0;ii<num_threads;ii++) t_dp_max[ii] =0;
+ for (int ii=0;ii<num_threads;ii++) t_sum_p[ii]=0;
+ for (int ii=0;ii<num_threads;ii++) t_sum_dp[ii]=0;
+#endif
+ }
+ #ifdef GS_OPEN_MP
+ #pragma omp parallel for schedule(dynamic,2000)
+ #endif
+
+ for (int ii=1; ii<=ncols; ii++) {
+ double** Acols = &(fullAcolumns[ii][0]); double* Avals = &(fullAvalues[ii][0]);
+ double dp = (((gsB[ii]-gsdV[ii]+Avals[0]*(*Acols[0])
+ +Avals[1]*(*Acols[1])
+ +Avals[2]*(*Acols[2])
+ +Avals[3]*(*Acols[3])) * Avals[4])
+ - gsP[ii])*relax;
+
+ gsP[ii]=dp+gsP[ii];
+ if (j2==0) {
+#ifdef GS_OPEN_MP
+ const int tn=omp_get_thread_num();
+ t_sum_dp[tn] += std::abs(dp);
+ t_dp_max[tn]=max(t_dp_max[tn], std::abs(dp));
+ t_p_max[tn]= max(t_p_max[tn], gsP[ii]);
+ t_sum_p[tn]+= std::abs(gsP[ii]);
+#else
+ dp_max = max(dp_max, std::abs(dp));
+ p_max = max(p_max, std::abs(gsP[ii]));
+ sum_p += std::abs(gsP[ii]);
+ sum_dp += std::abs(dp);
+#endif
+ }
+ }
+#ifdef GS_OPEN_MP
+ if (j2==0) {
+ for (int jj=0;jj<num_threads;jj++) p_max =max(p_max, t_p_max[jj]);
+ for (int jj=0;jj<num_threads;jj++) dp_max =max(dp_max, t_dp_max[jj]);
+ for (int jj=0;jj<num_threads;jj++) sum_p+=t_sum_p[jj];
+ for (int jj=0;jj<num_threads;jj++) sum_dp+=t_sum_dp[jj];
+ }
+#endif
+ if (j2==0) {
+ p_moy = sum_p/ncols;
+ dp_moy = sum_dp/ncols;
+ if (DEBUG_OUT) cerr <<"GS : j="<<j<<" p_moy="<<p_moy<<" dp_moy="<<dp_moy<<endl;
+ }
+#ifdef GS_OPEN_MP
+#pragma omp master
+#endif
+ j++;
+ } while ((dp_max/p_max) > tolerance && j<20000 /*&& ( dp_max > tolerance )*//* &&*/ /*( j<50 )*/);
+ CopyGsToCells();
+ if (j>=20000) cerr<<"GS did not converge in 20k iterations (maybe because the reference pressure is 0?)"<<endl;
+ if (DEBUG_OUT) cerr <<"GS iterations : "<<j-1<<endl;
+}
+
+template<class FlowType>
+void FlowBoundingSphereLinSolv<FlowType>::sort_v(int k1, int k2, int* is, double* ds){
+ for (int k=k1; k<k2; k++) {
+ int kk=k;
+ while (kk>=k1 && is[kk]>is[kk+1]) {
+ swap_fwd(is,kk);
+ swap_fwd(ds,kk);
+ --kk;}
+ }
+}
+
+template<class FlowType>
+int FlowBoundingSphereLinSolv<FlowType>::eigenSolve(Real dt)
+{
+#ifdef EIGENSPARSE_LIB
+ if (!isLinearSystemSet || (isLinearSystemSet && reApplyBoundaryConditions())) ncols = SetLinearSystem(dt);
+ CopyCellsToLin(dt);
+ //FIXME: we introduce new Eigen vectors, then we have to copy from/to c-arrays, can be optimized later
+ Eigen::VectorXd eb(ncols); Eigen::VectorXd ex(ncols);
+ for (int k=0; k<ncols; k++) eb[k]=T_bv[k];
+ if (!factorizedEigenSolver) {
+ eSolver.setMode(Eigen::CholmodSupernodalLLt);
+ openblas_set_num_threads(numFactorizeThreads);
+ eSolver.compute(A);
+ //Check result
+ if (eSolver.cholmod().status>0) {
+ cerr << "something went wrong in Cholesky factorization, use LDLt as fallback this time" << endl;
+ eSolver.setMode(Eigen::CholmodLDLt);
+ eSolver.compute(A);
+ }
+ factorizedEigenSolver = true;
+ }
+ openblas_set_num_threads(numSolveThreads);
+ ex = eSolver.solve(eb);
+ for (int k=0; k<ncols; k++) T_x[k]=ex[k];
+ CopyLinToCells();
+#else
+ cerr<<"Flow engine not compiled with eigen, nothing computed if useSolver=3"<<endl;
+#endif
+ return 0;
+}
+
+
+template<class FlowType>
+int FlowBoundingSphereLinSolv<FlowType>::TaucsSolve(Real dt)
+{
+#ifdef TAUCS_LIB
+ if (DEBUG_OUT) cerr <<endl<<"TAUCS solve"<<endl;
+ double t = taucs_ctime();//timer
+ double t2 = taucs_ctime();//global timer
+ if (!isLinearSystemSet || (isLinearSystemSet && reApplyBoundaryConditions())) {
+ ncols = SetLinearSystem(dt);
+ if (DEBUG_OUT) cerr << "Assembling the matrix : " << taucs_ctime()-t << endl; t = taucs_ctime();}
+
+ CopyCellsToLin(dt);
+ if (DEBUG_OUT) cerr << "Updating dv's (Yade->LinSolver) : " << taucs_ctime()-t << endl; t = taucs_ctime();
+ //taucs_logfile("stdout");//! VERY USEFULL!!!!!!!!!!! (disable to exclude output time from taucs_ctime() measurments)
+
+ taucs_double* x = &T_x[0];// the unknown vector to solve Ax=b
+ taucs_double* bod = &bodv[0];
+ taucs_double* xod = &xodv[0];
+
+ if (Fccs==NULL) {
+ if (DEBUG_OUT) cerr << "_entering taucs_" << endl;
+ // 1) Reordering
+ taucs_ccs_order(T_A, &perm, &invperm, (char*)"metis");
+ if (DEBUG_OUT) cerr << "_entering taucs_" << endl;
+ taucs_ccs_matrix* Aod;
+ if (DEBUG_OUT) cerr << "_entering taucs_" << endl;
+ Aod = taucs_ccs_permute_symmetrically(T_A, perm, invperm);
+ if (DEBUG_OUT) cerr << "Reordering : " << taucs_ctime()-t << endl; t = taucs_ctime();
+
+ // 2) Factoring
+ F = taucs_ccs_factor_llt_mf(Aod);
+ if (F==NULL) cerr<<"factorization failed"<<endl;
+ taucs_dccs_free(Aod); Aod=NULL;
+ //convert F to ccs for faster solve
+ Fccs = taucs_supernodal_factor_to_ccs(F);
+ //... then delete F
+ taucs_supernodal_factor_free(F); F=NULL;
+ if (DEBUG_OUT) cerr << "Factoring : " << taucs_ctime()-t << endl; t = taucs_ctime();
+ }
+ taucs_vec_permute(ncols, TAUCS_DOUBLE, &T_bv[0], bod, perm);
+ // 3) Back substitution and reodering the solution back
+ taucs_ccs_solve_llt(Fccs, xod, bod);//the ccs format (faster)
+// taucs_supernodal_solve_llt(F, xod, bod);//the blackbox format (slower)
+ if (DEBUG_OUT) cerr << "Solving : " << taucs_ctime()-t << endl; t = taucs_ctime();
+ t = taucs_ctime();
+ taucs_vec_ipermute(ncols, TAUCS_DOUBLE, xod, x, perm);
+// cerr << "Deordering : " << taucs_ctime()-t << endl; t = taucs_ctime();
+ // 4) Copy back to the triangulation
+ CopyLinToCells();
+// cerr << "Updating P (LinSolver->Yade) : " << taucs_ctime()-t << endl;
+ if (DEBUG_OUT) cerr << "Total TAUCS time ................ : " << taucs_ctime()-t2 << endl;
+#else
+ cerr<<"Flow engine not compiled with taucs, nothing computed if useSolver=1"<<endl;
+#endif
+ return 0;
+}
+template<class FlowType>
+int FlowBoundingSphereLinSolv<FlowType>::PardisoSolve(Real dt)
+{
+ cerr <<endl<<"PardisoSolve solve"<<endl;
+ #ifndef PARDISO
+ return 0;
+ #else
+ double iniT = taucs_ctime();
+
+ if (DEBUG_OUT) cerr << "_entering pardiso_" << endl;
+ /* Matrix data. */
+ double t = taucs_ctime();//timer
+ bool wasLSystemSet= isLinearSystemSet;
+ if (!isLinearSystemSet || (isLinearSystemSet && reApplyBoundaryConditions())) {
+ ncols = SetLinearSystem(dt);
+ if (DEBUG_OUT) cerr << "Assembling the matrix : " << taucs_ctime()-t << endl; t = taucs_ctime();}
+
+ if (DEBUG_OUT) cerr<<taucs_ctime()-t<<"s : set system"<<endl;
+ t=taucs_ctime();
+ ia = T_A->colptr;
+ ja = T_A->rowind;
+ a = T_A->values.d;
+ if (DEBUG_OUT) cerr<<taucs_ctime()-t<<"s : set system"<<endl;
+ if (!wasLSystemSet) for (int k=0; k<ncols; k++) sort_v(ia[k],ia[k+1]-1,ja,a);
+ if (DEBUG_OUT) cout<<taucs_ctime()-t<<"s for ordering CCS format"<<endl;
+ t=taucs_ctime();
+
+ nnz = ia[ncols];
+// int mtype = -2; /* Real symmetric matrix */
+ mtype = 2; /* Real symmetric positive def. matrix */
+ /* RHS and solution vectors. */
+ CopyCellsToLin(dt);
+ b = &T_bv[0];
+// P_x.resize(n);
+ x = &T_x[0];// the unknown vector to solve Ax=b
+ nrhs = 1; /* Number of right hand sides. */
+
+ error = 0;
+ solver = 0; /* use sparse direct solver */
+
+ /* Numbers of processors, value of OMP_NUM_THREADS */
+
+ if (!pardisoInitialized) {
+
+ var = getenv("OMP_NUM_THREADS");
+ if (var != NULL)
+ sscanf(var, "%d", &num_procs);
+ else {
+ num_procs=1;
+ cerr<<"Set environment OMP_NUM_THREADS to something. Pardiso needs it defined!"<<endl;
+ }
+ if (DEBUG_OUT) cerr<<taucs_ctime()-t<<"pardisoinit"<<endl;
+ F77_FUNC(pardisoinit)(pt, &mtype, &solver, iparm, dparm, &error);
+ if (DEBUG_OUT) cerr<<taucs_ctime()-t<<"pardisoinit'ed"<<endl;
+ pardisoInitialized=true;
+ if (error != 0) {
+ if (error == -10) printf("No license file found \n");
+ if (error == -11) printf("License is expired \n");
+ if (error == -12) printf("Wrong username or hostname \n");
+ return 1;}
+ iparm[2] = num_procs;
+ maxfct = 1; /* Maximum number of numerical factorizations. */
+ mnum = 1; /* Which factorization to use. */
+ msglvl = 0; /* Print statistical information */
+ error = 0; /* Initialize error flag */
+
+ /* .. Convert matrix from 0-based C-notation to Fortran 1-based */
+ /* notation. */
+ if (DEBUG_OUT) cout<<taucs_ctime()-t<<"tuning"<<endl;
+ t=taucs_ctime();
+ for (i = 0; i < ncols+1; i++) {
+ ia[i] += 1;
+ }
+ for (i = 0; i < nnz; i++) {
+ ja[i] += 1;
+ }
+ if (DEBUG_OUT) cout<<taucs_ctime()-t<<"s : Convert matrix from 0-based"<<endl;
+ t=taucs_ctime();
+ /* .. Reordering and Symbolic Factorization. This step also allocates */
+ /* all memory that is necessary for the factorization. */
+ phase = 11;
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &ncols, a, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, &ddum, &ddum, &error, dparm);
+ if (error != 0) {
+ printf("\nERROR during symbolic factorization: %d", error);
+ exit(1);
+ }
+ if (DEBUG_OUT) cout<<taucs_ctime()-t<<"s : Reordering and Symbolic Factorization"<<endl;
+ t=taucs_ctime();
+
+ /* .. Numerical factorization. */
+ phase = 22;
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &ncols, a, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, &ddum, &ddum, &error, dparm);
+
+ if (error != 0) {
+ printf("\nERROR during numerical factorization: %d", error);
+ exit(2);
+ }
+ if (DEBUG_OUT) cerr<<taucs_ctime()-t<<"s : Numerical factorization. "<<endl;
+ t=taucs_ctime();
+ }
+ /* .. Back substitution and iterative refinement. */
+ phase = 33;
+ iparm[7] = 0; /* Max numbers of iterative refinement steps. */
+ t=taucs_ctime();
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &ncols, a, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, b, x, &error, dparm);
+ if (error != 0) {
+ printf("\nERROR during solution: %d", error);
+ exit(3);
+ }
+ if (DEBUG_OUT) cerr<<taucs_ctime()-t<<"s : Back substitution and iterative refinement."<<endl;
+ t=taucs_ctime();
+ CopyLinToCells();
+ if (DEBUG_OUT) cerr<<taucs_ctime()-t<<"s : Copy back."<<endl;
+
+ if (wasLSystemSet){
+ pTime1N++; pTime1+=(taucs_ctime()-iniT);
+ } else {
+ pTime2N++; pTime2+=(taucs_ctime()-iniT);
+ }
+ if (pTimeInt>99) {
+ cout <<"Pardiso..... "<<pTime1/(double) pTime1N << " s/iter for "<<pTime1N<<"/"<<(pTime2N+pTime1N)<<" std iter., "<< pTime2/pTime2N <<" s/iter for "<<pTime2N<<"/"<<(pTime2N+pTime1N)<<" retriangulation iter."<<endl;
+ pTime1=0;pTime2=0;pTime1N=0;pTime2N=0;pTimeInt=0;}
+ pTimeInt++;
+ /* .. Convert matrix back to 0-based C-notation. */
+// for (i = 0; i < n+1; i++) {
+// ia[i] -= 1;
+// }
+// for (i = 0; i < nnz; i++) {
+// ja[i] -= 1;
+// }
+ return 0;
+ #endif
+}
+
+
+template<class FlowType>
+int FlowBoundingSphereLinSolv<FlowType>::PardisoSolveTest()
+{
+ #ifndef PARDISO
+ return 0;
+ #else
+ /* Matrix data. */
+ double t = taucs_ctime();//timer
+ bool wasLSystemSet= isLinearSystemSet;
+ int n = SetLinearSystem();
+// ncols=n;//for VectorizesGS
+ cout<<taucs_ctime()-t<<"s : set system"<<endl;
+ t=taucs_ctime();
+ int* ia = T_A->colptr;
+ int* ja = T_A->rowind;
+ double* a = T_A->values.d;
+
+ if (!wasLSystemSet) for (int k=0; k<n; k++) sort_v(ia[k],ia[k+1]-1,ja,a);
+ cout<<taucs_ctime()-t<<"s for ordering CCS format"<<endl;
+ t=taucs_ctime();
+
+ int nnz = ia[n];
+// int mtype = -2; /* Real symmetric matrix */
+ int mtype = 2; /* Real symmetric positive def. matrix */
+ /* RHS and solution vectors. */
+ double* b = &T_b[0];
+ P_x.resize(n);
+ double* x = &P_x[0];// the unknown vector to solve Ax=b
+ int nrhs = 1; /* Number of right hand sides. */
+
+ /* Internal solver memory pointer pt, */
+ /* 32-bit: int pt[64]; 64-bit: long int pt[64] */
+ /* or void *pt[64] should be OK on both architectures */
+ void *pt[64];
+ /* Pardiso control parameters. */
+ int iparm[64];
+ double dparm[64];
+ int maxfct, mnum, phase, error, msglvl, solver;
+ /* Number of processors. */
+ int num_procs;
+ /* Auxiliary variables. */
+ char *var;
+ int i;
+
+ double ddum; /* Double dummy */
+ int idum; /* Integer dummy. */
+ /* .. Setup Pardiso control parameters. */
+ error = 0;
+ solver = 0; /* use sparse direct solver */
+
+ /* Numbers of processors, value of OMP_NUM_THREADS */
+ var = getenv("OMP_NUM_THREADS");
+ if (var != NULL)
+ sscanf(var, "%d", &num_procs);
+ else {
+ num_procs=1;
+ printf("Set environment OMP_NUM_THREADS to something. Pardiso needs it defined. \n");
+// exit(0);
+ }
+
+ F77_FUNC(pardisoinit)(pt, &mtype, &solver, iparm, dparm, &error);
+ if (error != 0)
+ {
+ if (error == -10)
+ printf("No license file found \n");
+ if (error == -11)
+ printf("License is expired \n");
+ if (error == -12)
+ printf("Wrong username or hostname \n");
+ return 1;
+ }
+ iparm[2] = num_procs;
+ maxfct = 1; /* Maximum number of numerical factorizations. */
+ mnum = 1; /* Which factorization to use. */
+ msglvl = 0; /* Print statistical information */
+ error = 0; /* Initialize error flag */
+
+ /* .. Convert matrix from 0-based C-notation to Fortran 1-based */
+ /* notation. */
+ cout<<taucs_ctime()-t<<"tuning"<<endl;
+ t=taucs_ctime();
+ for (i = 0; i < n+1; i++) {
+ ia[i] += 1;
+ }
+ for (i = 0; i < nnz; i++) {
+ ja[i] += 1;
+ }
+ cout<<taucs_ctime()-t<<"s : Convert matrix from 0-based"<<endl;
+ t=taucs_ctime();
+ /* .. Reordering and Symbolic Factorization. This step also allocates */
+ /* all memory that is necessary for the factorization. */
+ phase = 11;
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &n, a, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, &ddum, &ddum, &error, dparm);
+ if (error != 0) {
+ printf("\nERROR during symbolic factorization: %d", error);
+ exit(1);
+ }
+ cout<<taucs_ctime()-t<<"s : Reordering and Symbolic Factorization"<<endl;
+ t=taucs_ctime();
+
+ /* .. Numerical factorization. */
+ phase = 22;
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &n, a, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, &ddum, &ddum, &error, dparm);
+
+ if (error != 0) {
+ printf("\nERROR during numerical factorization: %d", error);
+ exit(2);
+ }
+ cerr<<taucs_ctime()-t<<"s : Numerical factorization. "<<endl;
+ t=taucs_ctime();
+
+ /* .. Back substitution and iterative refinement. */
+ phase = 33;
+ iparm[7] = 0; /* Max numbers of iterative refinement steps. */
+ t=taucs_ctime();
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &n, a, ia, ja, &idum, &nrhs,
+ iparm, &msglvl, b, x, &error, dparm);
+ if (error != 0) {
+ printf("\nERROR during solution: %d", error);
+ exit(3);
+ }
+ cerr<<taucs_ctime()-t<<"s : Back substitution and iterative refinement."<<endl;
+ t=taucs_ctime();
+
+ /* .. Convert matrix back to 0-based C-notation. */
+ for (i = 0; i < n+1; i++) {
+ ia[i] -= 1;
+ }
+ for (i = 0; i < nnz; i++) {
+ ja[i] -= 1;
+ }
+ /* .. Termination and release of memory. */
+ phase = -1; /* Release internal memory. */
+// F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+// &n, &ddum, ia, ja, &idum, &nrhs,
+// iparm, &msglvl, &ddum, &ddum, &error, dparm);
+ //Release only what has to be
+ F77_FUNC(pardiso)(pt, &maxfct, &mnum, &mtype, &phase,
+ &n, &ddum, NULL, NULL, &idum, &nrhs,
+ iparm, &msglvl, &ddum, &ddum, &error, dparm);
+ return 0;
+ #endif
+}
+template<class FlowType>
+int FlowBoundingSphereLinSolv<FlowType>::TaucsSolveTest()
+{
+#ifdef TAUCS_LIB
+ cout <<endl<<"TAUCS solve test"<<endl;
+
+ double t = taucs_ctime();//timer
+ ncols = SetLinearSystem();
+
+//taucs_logfile("stdout");//! VERY USEFULL!!!!!!!!!!! (disable to exclude output time from taucs_ctime() measurments)
+
+// allocate TAUCS solution vector
+ T_x.resize(ncols);
+ double* x = &*T_x.begin();// the unknown vector to solve Ax=b
+ cout << "Assembling the matrix2 : " << taucs_ctime()-t << endl;
+ t =taucs_ctime();
+// solve the linear system
+ void* F = NULL;
+
+//Allocate reoredered x and b
+ vector <double> bodv(ncols);
+ taucs_double* bod = &*bodv.begin();
+ vector <double> xodv(ncols);
+ taucs_double* xod = &*xodv.begin();
+
+
+ int* perm;
+ int* invperm;
+ taucs_ccs_matrix* Aod;
+
+ t = taucs_ctime();
+ double t2 = taucs_ctime();
+// 1) Reordering
+ taucs_ccs_order(T_A, &perm, &invperm, "metis");
+ Aod = taucs_ccs_permute_symmetrically(T_A, perm, invperm);
+ taucs_vec_permute(ncols, TAUCS_DOUBLE, &T_b[0], bod, perm);
+ cout << "Reordering : " << taucs_ctime()-t << endl;
+ t = taucs_ctime();
+// 2) Factoring
+
+ F = taucs_ccs_factor_llt_mf(Aod);
+//F = taucs_ccs_factor_llt_mf(T_A);
+ cout << "Factoring : " << taucs_ctime()-t << endl;
+ t = taucs_ctime();
+// 3) Back substitution and reodering the solution back
+double t4 = taucs_ctime();
+// for (int k=0;k<10;k++){
+
+ taucs_supernodal_solve_llt(F, xod, bod);
+// cout << "B3) Solving : " << taucs_ctime()-t << endl;
+ t = taucs_ctime();
+ taucs_vec_ipermute(ncols, TAUCS_DOUBLE, xod, x, perm);
+// cout << "B4) Deordering : " << taucs_ctime()-t << endl;
+ t = taucs_ctime();
+// }
+ double T4=taucs_ctime()-t4;
+ cout << "Solving : " << T4 << endl;
+ cout << "Low level reordered total time : " << taucs_ctime()-t2 << endl;
+ t2 = taucs_ctime();
+ taucs_supernodal_factor_free(F);
+ taucs_ccs_free(Aod);
+
+/// Using TAUCS inverse factoring
+// t = taucs_ctime(); t2 = taucs_ctime();
+// F=taucs_ccs_factor_xxt(Aod);
+// cout << "C1_) inverse factoring (reordered) : " << taucs_ctime()-t << endl; t = taucs_ctime();
+// taucs_ccs_solve_xxt (F,xod,bod);///REALLY TOO SLOW!!
+// cout << "C2) solve : " << taucs_ctime()-t << endl; t = taucs_ctime();
+// taucs_vec_ipermute(ncols, TAUCS_DOUBLE, xod, x, perm);
+// cout << "C3) reordering : " << taucs_ctime()-t << endl; t = taucs_ctime();
+// cout << "Inverse factoring total time : " << taucs_ctime()-t2 << endl;
+// taucs_linsolve(NULL, &F, 0, NULL, NULL, NULL, NULL);
+// taucs_ccs_free(Aod);
+
+// ofstream file("result.dat");
+// if (!file.is_open()) cout << "problem opening file";
+// for (unsigned i = 0; i < M ; i++) {
+// for (unsigned j = 0; j < L ; j++) file << i << " "<< j << " "<< x[index(i,j)] << endl;
+// file << endl << endl;
+// }
+
+
+// const Finite_cells_iterator cell_end = T[currentTes].Triangulation().finite_cells_end();
+// if (DEBUG_OUT) for (Finite_cells_iterator cell = T[currentTes].Triangulation().finite_cells_begin(); cell != cell_end; cell++)
+// {
+//
+// if (cell->info().index>0) {
+// cerr <<"ix "<< cell->info().index<<" ";
+// cerr <<cell->info().p()<<" "<<gsP[cell->info().index]<<" "<<P_x[max((unsigned int)0,cell->info().index-1)] <<" "<<x[max((unsigned int)0,cell->info().index-1)]<<" err="<<100*(cell->info().p()-x[max((unsigned int)0,cell->info().index-1)])/x[max((unsigned int)0,cell->info().index-1)]<<"%"<<endl;
+// }
+// else if (cell->info().index==0) {
+// cerr <<"ix "<< cell->info().index<<" ";
+// cerr <<cell->info().p()<<" "<<"NaN"<<endl;
+// }
+// }
+
+#endif
+ return 0;
+}
+
+
+} //namespace CGT
+
+#endif //FLOW_ENGINE
\ No newline at end of file
=== added file 'lib/triangulation/PeriodicFlowLinSolv.hpp'
--- lib/triangulation/PeriodicFlowLinSolv.hpp 1970-01-01 00:00:00 +0000
+++ lib/triangulation/PeriodicFlowLinSolv.hpp 2013-06-27 11:29:17 +0000
@@ -0,0 +1,38 @@
+/*************************************************************************
+* Copyright (C) 2010 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. *
+*************************************************************************/
+
+#ifndef _PERIODICFLOWLINSOLV_H
+#define _PERIODICFLOWLINSOLV_H
+#include "FlowBoundingSphere.hpp"
+
+#ifdef FLOW_ENGINE
+
+using namespace std;
+
+namespace CGT {
+
+typedef FlowBoundingSphereLinSolv<PeriodicFlow> LinSolver;
+
+class PeriodicFlowLinSolv : public LinSolver
+{
+public:
+ typedef PeriodicFlow FlowType;
+ vector<int> indices;//redirection vector containing the rank of cell so that T_cells[indices[cell->info().index]]=cell
+
+ virtual ~PeriodicFlowLinSolv();
+ PeriodicFlowLinSolv();
+
+ ///Linear system solve
+ virtual int SetLinearSystem(Real dt=0);
+ virtual int SetLinearSystemFullGS(Real dt=0);
+};
+
+} //namespace CGTF
+
+#endif //FLOW_ENGINE
+
+#endif
=== added file 'lib/triangulation/PeriodicFlowLinSolv.ipp'
--- lib/triangulation/PeriodicFlowLinSolv.ipp 1970-01-01 00:00:00 +0000
+++ lib/triangulation/PeriodicFlowLinSolv.ipp 2013-06-27 11:29:17 +0000
@@ -0,0 +1,325 @@
+/*************************************************************************
+* Copyright (C) 2010 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 FLOW_ENGINE
+
+
+#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#include <CGAL/Width_3.h>
+#include <iostream>
+#include <fstream>
+#include <new>
+#include <utility>
+#include "vector"
+#include <assert.h>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <omp.h>
+
+
+namespace CGT
+{
+
+using namespace std;
+
+#ifdef PARDISO
+#ifdef AIX
+#define F77_FUNC(func) func
+#else
+#define F77_FUNC(func) func ## _
+#endif
+/* PARDISO prototype. */
+extern "C" int F77_FUNC(pardisoinit)
+ (void *, int *, int *, int *, double *, int *);
+
+extern "C" int F77_FUNC(pardiso)
+ (void *, int *, int *, int *, int *, int *,
+ double *, int *, int *, int *, int *, int *,
+ int *, double *, double *, int *, double *);
+#endif
+
+PeriodicFlowLinSolv::~PeriodicFlowLinSolv()
+{
+}
+
+PeriodicFlowLinSolv::PeriodicFlowLinSolv(): LinSolver() {}
+
+int PeriodicFlowLinSolv::SetLinearSystem(Real dt)
+{
+//WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition of
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ vector<int> clen;
+ vector<int> is;
+ vector<int> js;
+ vector<double> vs;
+ if (!areCellsOrdered) {
+ T_nnz=0;
+ ncols=0;
+ T_cells.clear();
+ T_index=0;
+ unsigned int maxindex = 0;
+ //FIXME: this is way too large since many cells will be ghosts
+ T_cells.resize(Tri.number_of_finite_cells()+1);
+// indices.resize(Tri.number_of_finite_cells()+1);
+ ///Ordered cells
+ orderedCells.clear();
+ const Finite_cells_iterator cell_end = Tri.finite_cells_end();
+ for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ if (cell->info().Pcondition || cell->info().isGhost) continue;
+ orderedCells.push_back(cell);
+// T_cells[++ncols]=cell;
+// indices[cell->info().index]=ncols;
+ ++ncols;
+ T_cells[cell->info().index]=cell;
+ maxindex=max(maxindex,cell->info().index);
+ }
+// spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort());//FIXME: ordering will not work with the new "indices", which needs reordering to
+ T_cells.resize(ncols+1);
+// indices.resize(maxindex+1);
+ isLinearSystemSet=false;
+ areCellsOrdered=true;
+ }
+ if (!isLinearSystemSet) {
+ int n = 3*(ncols+1);//number of non-zero in triangular matrix
+ is.resize(n);
+ js.resize(n);
+ vs.resize(n);
+ T_x.resize(ncols);
+ T_b.resize(ncols);
+ T_bv.resize(ncols);
+ bodv.resize(ncols);
+ xodv.resize(ncols);
+// gsB.resize(ncols+1);
+ T_cells.resize(ncols+1);
+ T_nnz=0;}
+ for (int kk=0; kk<ncols;kk++) T_b[kk]=0;
+ ///Ordered cells
+ int nIndex=0; Cell_handle neighbour_cell;
+ for (int i=0; i<ncols; i++)
+ {
+ Finite_cells_iterator& cell = orderedCells[i];
+ ///Non-ordered cells
+// for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+// if (!cell->info().Pcondition && !cell->info().isGhost)
+ {
+ const int& index=cell->info().index;
+ const Cell_Info& cInfo = cell->info();
+ if (!isLinearSystemSet) {
+ //Add diagonal term
+ is[T_nnz] = index;
+ js[T_nnz] = index;
+ vs[T_nnz] = (cInfo.k_norm())[0]+ (cInfo.k_norm())[1]+ (cInfo.k_norm())[2]+ (cInfo.k_norm())[3];
+ if (vs[T_nnz]<0) cerr<<"!!!! WTF !!!"<<endl;
+ if (fluidBulkModulus>0) vs[T_nnz] += (1.f/(dt*fluidBulkModulus*cInfo.invVoidVolume()));
+ T_nnz++;
+ }
+ for (int j=0; j<4; j++) {
+ neighbour_cell = cell->neighbor(j);
+ if (Tri.is_infinite(neighbour_cell)) continue;
+ Cell_Info& nInfo = neighbour_cell->info();
+ nIndex=nInfo.index;
+ if (nIndex==index) {
+ cerr<<"ERROR: nIndex==index, the cell is neighbour to itself"<< index<<endl;
+ errorCode=3;}
+
+ if (nInfo.Pcondition) T_b[index-1]+=cInfo.k_norm()[j]*nInfo.shiftedP();
+ else {
+ if (!isLinearSystemSet && index>nIndex) {
+ is[T_nnz] = index;
+ js[T_nnz] = nIndex;
+ vs[T_nnz] = - (cInfo.k_norm())[j];
+ if (vs[T_nnz]>0) cerr<<"!!!! WTF2 !!!"<<endl;
+ T_nnz++;}
+ if (nInfo.isGhost) T_b[index-1]+=cInfo.k_norm()[j]*nInfo.pShift();
+ }
+ }
+ }
+ }
+ if (!isLinearSystemSet) {
+ if (useSolver==1 || useSolver==2){
+ #ifdef TAUCS_LIB
+ clen.resize(ncols+1);
+ T_jn.resize(ncols+1);
+ T_A->colptr = &T_jn[0];
+ T_ia.resize(T_nnz);
+ T_A->rowind = &T_ia[0];
+ T_A->flags = (TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
+ T_an.resize(T_nnz);
+ T_A->values.d = &T_an[0];
+ T_A->n = ncols;
+ T_A->m = ncols;
+ int i,j,k;
+ for (j=0; j<ncols; j++) clen[j] = 0;
+ for (k=0; k<T_nnz; k++) {
+ i = is[k]-1; /* make it 1-based */
+ j = js[k]-1; /* make it 1-based */
+ (clen[j])++;
+ }
+ /* now compute column pointers */
+ k = 0;
+ for (j=0; j<ncols; j++) {
+ int tmp;
+ tmp = clen[j];
+ clen[j] = (T_A->colptr[j]) = k;
+ k += tmp;
+ }
+ clen[ncols] = (T_A->colptr[ncols]) = k;
+
+ /* now read matrix into data structure */
+ for (k=0; k<T_nnz; k++) {
+ i = is[k] - 1; /* make it 1-based */
+ j = js[k] - 1; /* make it 1-based */
+ assert(i < ncols);
+ assert(j < ncols);
+ (T_A->taucs_values)[clen[j]] = vs[k];
+ (T_A->rowind)[clen[j]] = i;
+ clen[j] ++;
+ }
+ #else
+ cerr<<"yade compiled without Taucs, FlowEngine.useSolver="<< useSolver <<" not supported"<<endl;
+ #endif //TAUCS_LIB
+ } else if (useSolver==3){
+ #ifdef EIGENSPARSE_LIB
+// //here the matrix can be exported in in MatrixMarket format for experiments
+// static int mn=0;
+// ofstream file; ofstream file2; ofstream file3;
+// stringstream ss,ss2,ss3;
+// ss<<"matrix"<<mn;
+// ss3<<"matrixf2"<<mn;
+// ss2<<"matrixf"<<mn++;
+// file.open(ss.str().c_str());
+// file2.open(ss2.str().c_str());
+// file3.open(ss3.str().c_str());
+// file <<"%%MatrixMarket matrix coordinate real symmetric"<<endl;
+// file2 <<"%%MatrixMarket matrix coordinate real symmetric"<<endl;
+// file3 <<"%%MatrixMarket matrix coordinate real symmetric"<<endl;
+// file <<ncols<<" "<<ncols<<" "<<T_nnz<<" -1"<<endl;
+// file2 <<ncols<<" "<<ncols<<" "<<T_nnz<<" -1"<<endl;
+// file3 <<ncols<<" "<<ncols<<" "<<T_nnz<<" -1"<<endl;
+ tripletList.clear(); tripletList.resize(T_nnz);
+ for(int k=0;k<T_nnz;k++) {
+ tripletList[k]=ETriplet(is[k]-1,js[k]-1,vs[k]);
+// file<<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
+// if (is[k]==js[k]) file2<<is[k]-1<<" "<<js[k]-1<<" "<<1.0001*vs[k]<<endl;
+// else file2<<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
+// if (is[k]==js[k]) file3<<is[k]-1<<" "<<js[k]-1<<" "<<1.00000001*vs[k]<<endl;
+// else file3<<is[k]-1<<" "<<js[k]-1<<" "<<vs[k]<<endl;
+ }
+ A.resize(ncols,ncols);
+ A.setFromTriplets(tripletList.begin(), tripletList.end());
+// file << A;
+// file.close();
+ #else
+ cerr<<"yade compiled without CHOLMOD, FlowEngine.useSolver="<< useSolver <<" not supported"<<endl;
+ #endif
+ }
+ isLinearSystemSet=true;
+ }
+ return ncols;
+}
+
+
+
+/// For Gauss Seidel, we need the full matrix
+
+int PeriodicFlowLinSolv::SetLinearSystemFullGS(Real dt)
+{
+ //WARNING : boundary conditions (Pcondition, p values) must have been set for a correct definition
+ RTriangulation& Tri = T[currentTes].Triangulation();
+ int n_cells=Tri.number_of_finite_cells();
+
+ if (!isFullLinearSystemGSSet){
+ T_cells.clear();
+ T_index=0;//FIXME : no need to clear if we don't re-triangulate
+ T_nnz=0;
+ ncols=0;
+ const Finite_cells_iterator cell_end = Tri.finite_cells_end();
+ orderedCells.clear();
+ T_cells.resize(n_cells+1);
+ for (Finite_cells_iterator cell = Tri.finite_cells_begin(); cell != cell_end; cell++) {
+ if (cell->info().Pcondition || cell->info().isGhost) continue;
+ ++ncols;
+ T_cells[cell->info().index]=cell;
+ }
+// spatial_sort(orderedCells.begin(),orderedCells.end(), CellTraits_for_spatial_sort());
+ gsP.resize(ncols+1);
+ gsB.resize(ncols+1);
+ T_b.resize(ncols+1);
+ gsdV.resize(ncols+1);
+ fullAcolumns.resize(ncols+1);
+ fullAvalues.resize(ncols+1);
+ T_cells.resize(ncols+1);
+ for (int k=0; k<=ncols;k++) {
+ fullAcolumns[k].resize(4);
+ fullAvalues[k].resize(5);
+ gsdV[k]=0;
+ }
+ gsP[0]=0;
+ areCellsOrdered=true;
+ }
+ for (int k=0; k<=ncols;k++) gsB[k]=0;
+ if (!isFullLinearSystemGSSet){
+ ///Ordered cells
+ for (int i=1; i<=ncols; i++)
+ {
+ Cell_handle& cell = T_cells[i];
+ ///Non-ordered cells
+ if (!cell->info().Pcondition && !cell->info().isGhost) {
+ //Add diagonal term
+ fullAvalues[cell->info().index][4] = 1.f/((cell->info().k_norm())[0]+ (cell->info().k_norm())[1]+ (cell->info().k_norm())[2]+ (cell->info().k_norm())[3] + (fluidBulkModulus>0? 1.f/(dt*fluidBulkModulus*cell->info().invVoidVolume()):0));
+ //DUMP
+// cout<< cell->info().index<<" "<< cell->info().index<<" "<<fullAvalues[cell->info().index][4] <<endl;
+ T_nnz++;
+ for (int j=0; j<4; j++) {
+ Cell_handle neighbour_cell = cell->neighbor(j);
+ const Cell_Info& nInfo = neighbour_cell->info();
+ Cell_Info& cInfo = cell->info();
+ if (Tri.is_infinite(neighbour_cell)) {
+ fullAvalues[cInfo.index][j] = 0;
+ //We point to the pressure of the adjacent cell. If the cell is ghost, then it has the index of the real one, and then the pointer is correct
+ fullAcolumns[cInfo.index][j] = &gsP[0];
+ continue;}
+ if (!nInfo.Pcondition) {
+ ++T_nnz;
+ fullAvalues[cInfo.index][j] = (cInfo.k_norm())[j];
+ fullAcolumns[cInfo.index][j] = &gsP[nInfo.index];
+ //DUMP
+// cout<< cInfo.index<<" "<< nInfo.index<<" "<<fullAvalues[cInfo.index][j] <<endl;
+ //if the adjacent cell is ghost, we account for the pressure shift in the RHS
+ if (nInfo.isGhost){
+ gsB[cInfo.index]+=cInfo.k_norm()[j]*nInfo.pShift();
+ }
+ } else {
+ fullAvalues[cInfo.index][j] = 0;
+ fullAcolumns[cInfo.index][j] = &gsP[0];
+ gsB[cInfo.index]+=cInfo.k_norm()[j]*nInfo.shiftedP();
+ }
+ }
+ }
+ }
+ } else for (int i=1; i<=ncols; i++)
+ {
+ Cell_handle& cell = T_cells[i];
+ ///Non-ordered cells
+ if (!cell->info().Pcondition && !cell->info().isGhost) {
+ for (int j=0; j<4; j++) {
+ Cell_handle neighbour_cell = cell->neighbor(j);
+ const Cell_Info& nInfo = neighbour_cell->info();
+ Cell_Info& cInfo = cell->info();
+ if (!nInfo.Pcondition) {
+ if (nInfo.isGhost) gsB[cInfo.index]+=cInfo.k_norm()[j]*nInfo.pShift();
+ } else gsB[cInfo.index]+=cInfo.k_norm()[j]*nInfo.shiftedP();
+ }
+ }
+ }
+ isFullLinearSystemGSSet=true;
+ return ncols;
+}
+
+} //namespace CGT
+
+#endif //FLOW_ENGINE
\ No newline at end of file
=== modified file 'pkg/dem/InelastCohFrictMat.hpp'
--- pkg/dem/InelastCohFrictMat.hpp 2012-12-14 14:38:44 +0000
+++ pkg/dem/InelastCohFrictMat.hpp 2013-06-25 08:02:13 +0000
@@ -1,37 +1,50 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
#pragma once
-#include "CohFrictMat.hpp"
-
-
-class InelastCohFrictMat : public CohFrictMat
+#include "../common/ElastMat.hpp"
+
+class InelastCohFrictMat : public FrictMat
{
public :
virtual ~InelastCohFrictMat ();
/// Serialization
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictMat,CohFrictMat,"",
- ((Real,eTT,0.0,,"Tension elasticity modulus"))
- ((Real,eTC,0.0,,"Compresion elasticity modulus"))
- ((Real,eB,0.0,,"bending elasticity modulus"))
- ((Real,gTw,0.0,,"twist elasticity modulus"))
- ((Real,sigmaB,0.0,,"max. Bending stress"))
- ((Real,sigmaTw,0.0,,"max. Twist stress"))
- ((Real,creepT,0.0,,"Tension/compression creeping"))
- ((Real,creepB,0.0,,"bending creeping"))
- ((Real,creepTw,0.0,,"twist creeping"))
- ((Real,unloadT,0.0,,"Tension/compression plastic unload"))
- ((Real,unloadB,0.0,,"bending plastic unload"))
- ((Real,unloadTw,0.0,,"twist plastic unload"))
- ((Real,disElT,0.0,,"max elastic displacement on Traction"))
- ((Real,disElC,0.0,,"max elastic displacement on Compresion"))
- ((Real,epsilonMaxT,0.0,,"Maximal plastic strain Tension"))
- ((Real,epsilonMaxC,0.0,,"Maximal plastic strain compression"))
- ((Real,phiBMax,0.0,,"Maximal plastic bending strain"))
- ((Real,phiTwMax,0.0,,"Maximal plastic bending strain")),
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictMat,FrictMat,"",
+ ((Real,tensionModulus,0.0,,"Tension elasticity modulus"))
+ ((Real,compressionModulus,0.0,,"Compresion elasticity modulus"))
+ ((Real,shearModulus,0.0,,"shear elasticity modulus"))
+ ((Real,alphaKr,2.0,,"Dimensionless coefficient used for the rolling stiffness."))
+ ((Real,alphaKtw,2.0,,"Dimensionless coefficient used for the twist stiffness."))
+
+ ((Real,nuBending,0.0,,"Bending elastic stress limit"))
+ ((Real,nuTwist,0.0,,"Twist elastic stress limit"))
+ ((Real,sigmaTension,0.0,,"Tension elastic stress limit"))
+ ((Real,sigmaCompression,0.0,,"Compression elastic stress limit"))
+ ((Real,shearCohesion,0.0,,"Shear elastic stress limit"))
+
+ ((Real,creepTension,0.0,,"Tension/compression creeping coefficient. Usual values between 0 and 1."))
+ ((Real,creepBending,0.0,,"Bending creeping coefficient. Usual values between 0 and 1."))
+ ((Real,creepTwist,0.0,,"Twist creeping coefficient. Usual values between 0 and 1."))
+
+ ((Real,unloadTension,0.0,,"Tension/compression plastic unload coefficient. Usual values between 0 and +infinity."))
+ ((Real,unloadBending,0.0,,"Bending plastic unload coefficient. Usual values between 0 and +infinity."))
+ ((Real,unloadTwist,0.0,,"Twist plastic unload coefficient. Usual values between 0 and +infinity."))
+
+ ((Real,epsilonMaxTension,0.0,,"Maximal plastic strain tension"))
+ ((Real,epsilonMaxCompression,0.0,,"Maximal plastic strain compression"))
+ ((Real,etaMaxBending,0.0,,"Maximal plastic bending strain"))
+ ((Real,etaMaxTwist,0.0,,"Maximal plastic twist strain")),
createIndex();
);
/// Indexable
- REGISTER_CLASS_INDEX(InelastCohFrictMat,CohFrictMat);
+ REGISTER_CLASS_INDEX(InelastCohFrictMat,FrictMat);
};
REGISTER_SERIALIZABLE(InelastCohFrictMat);
=== modified file 'pkg/dem/InelastCohFrictPhys.cpp'
--- pkg/dem/InelastCohFrictPhys.cpp 2012-12-14 14:38:44 +0000
+++ pkg/dem/InelastCohFrictPhys.cpp 2013-06-25 08:02:13 +0000
@@ -1,13 +1,12 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
#include "InelastCohFrictPhys.hpp"
-
-
-void InelastCohFrictPhys::SetBreakingState()
-{
- cohesionBroken = true;
- normalAdhesion = 0;
- shearAdhesion = 0;
-}
-
InelastCohFrictPhys::~InelastCohFrictPhys()
{
}
=== modified file 'pkg/dem/InelastCohFrictPhys.hpp'
--- pkg/dem/InelastCohFrictPhys.hpp 2013-05-14 21:24:11 +0000
+++ pkg/dem/InelastCohFrictPhys.hpp 2013-06-25 08:02:13 +0000
@@ -1,53 +1,73 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
#pragma once
-#include "CohFrictPhys.hpp"
-
-
-class InelastCohFrictPhys : public CohFrictPhys
+#include "FrictPhys.hpp"
+
+
+class InelastCohFrictPhys : public FrictPhys
{
public :
virtual ~InelastCohFrictPhys();
- void SetBreakingState ();
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictPhys,CohFrictPhys,"",
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(InelastCohFrictPhys,FrictPhys,"",
((bool,cohesionBroken,false,,"is cohesion active? will be set false when a fragile contact is broken"))
- ((bool,fragile,false,,"do cohesion disapear when contact strength is exceeded?"))
- ((Real,knT,0,,"tension stiffness [N.m/rad]"))
- ((Real,knC,0,,"compression stiffness [N.m/rad]"))
- ((Real,kt,0,,"twist shear stiffness [N.m/rad]"))
- ((Real,ks,0,,"shear stiffness [N.m/rad]"))
- ((Real,kr,0,,"bending stiffness [N.m/rad]"))
- ((Real,maxElastB,0.0,,"max. elastic Bending, zero if non elasto-plastic behavour "))
- ((Real,maxElastTw,0.0,,"max. elastic Twist, zero if non elasto-plastic behavour "))
- ((Real,dElT,0.0,,"Max elastic displacement Tension"))
- ((Real,dElC,0.0,,"Max elastic displacement Compression"))
- ((Real,crpT,0.0,,"Bending creep"))
- ((Real,crpB,0.0,,"Tension crepp"))
- ((Real,crpTw,0.0,,"Twist creep "))
- ((Real,epsMaxT,0.0,,"Maximal plastic strain Tension"))
- ((Real,epsMaxC,0.0,,"Maximal plastic strain compression"))
- ((Real,phBMax,0.0,,"Maximal plastic bending strain"))
- ((Real,phTwMax,0.0,,"Maximal plastic bending strain"))
- ((Real,unldT,0.0,,"Tension/compression plastic unload"))
- ((Real,unldB,0.0,,"bending plastic unload"))
- ((Real,unldTw,0.0,,"twist plastic unload"))
- ((bool,isBrokenB,false,,"true if bend plastic fracture achieved"))
- ((bool,isBrokenT,false,,"true if Traction or compression plastic fracture achieved"))
- ((bool,isBrokenTw,false,,"true if twist plastic fracture achieved"))
- ((bool,unloadedT,false,,"true if unload in Traction plastic deformation"))
- ((bool,unloadedC,false,,"true if unload in compression plastic deformation"))
- ((bool,unloadedB,false,,"true if unload in bending plastic deformation"))
- ((bool,unloadedTw,false,,"true if unload in twist plastic deformation"))
- ((Real,normalAdhesion,0,,"tensile strength"))
- ((Real,shearAdhesion,0,,"cohesive part of the shear strength (a frictional term might be added depending on :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::always_use_moment_law`)"))
- ((Real,unp,0,,"plastic normal displacement, only used for tensile behaviour and if :yref:`CohFrictPhys::fragile`=false."))
- ((bool,onPlastic,false,,"true if plastic behaviour achieved, for creeping calculation"))
+
+ ((Real,knT,0,,"tension stiffness"))
+ ((Real,knC,0,,"compression stiffness"))
+ ((Real,ktw,0,,"twist shear stiffness"))
+ ((Real,ks,0,,"shear stiffness"))
+ ((Real,kr,0,,"bending stiffness"))
+
+ ((Real,maxElB,0.0,,"Maximum bending elastic moment."))
+ ((Real,maxElTw,0.0,,"Maximum twist elastic moment."))
+ ((Real,maxElT,0.0,,"Maximum tension elastic force."))
+ ((Real,maxElC,0.0,,"Maximum compression elastic force."))
+ ((Real,shearAdhesion,0,,"Maximum elastic shear force (cohesion)."))
+
+ ((Real,kTCrp,0.0,,"Tension/compression creep stiffness"))
+ ((Real,kRCrp,0.0,,"Bending creep stiffness"))
+ ((Real,kTwCrp,0.0,,"Twist creep stiffness"))
+
+ ((Real,kTUnld,0.0,,"Tension/compression plastic unload stiffness"))
+ ((Real,kRUnld,0.0,,"Bending plastic unload stiffness"))
+ ((Real,kTwUnld,0.0,,"Twist plastic unload stiffness"))
+
+ ((Real,maxExten,0.0,,"Plastic failure extension (stretching)."))
+ ((Real,maxContract,0.0,,"Plastic failure contraction (shrinkage)."))
+ ((Real,maxBendMom,0.0,,"Plastic failure bending moment."))
+ ((Real,maxTwist,0.0,,"Plastic failure twist angle"))
+
+ ((bool,isBroken,false,,"true if compression plastic fracture achieved"))
+
+ ((Real,unp,0,,"plastic normal penetration depth describing the equilibrium state."))
+ ((Real,twp,0,,"plastic twist penetration depth describing the equilibrium state."))
+
+ ((bool,onPlastB,false,Attr::readonly,"true if plasticity achieved on bending"))
+ ((bool,onPlastTw,false,Attr::readonly,"true if plasticity achieved on twisting"))
+ ((bool,onPlastT,false,Attr::readonly,"true if plasticity achieved on traction"))
+ ((bool,onPlastC,false,Attr::readonly,"true if plasticity achieved on compression"))
+
+ ((Vector2r,maxCrpRchdT,Vector2r(0,0),Attr::readonly,"maximal extension reached on plastic deformation. maxCrpRchdT[0] stores un and maxCrpRchdT[1] stores Fn."))
+ ((Vector2r,maxCrpRchdC,Vector2r(0,0),Attr::readonly,"maximal compression reached on plastic deformation. maxCrpRchdC[0] stores un and maxCrpRchdC[1] stores Fn."))
+ ((Vector2r,maxCrpRchdTw,Vector2r(0,0),Attr::readonly,"maximal twist reached on plastic deformation. maxCrpRchdTw[0] stores twist angle and maxCrpRchdTw[1] stores twist moment."))
+ ((Vector3r,maxCrpRchdB,Vector3r(0,0,0),Attr::readonly,"maximal bending moment reached on plastic deformation."))
+
+ ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::readonly),"Twist moment"))
+ ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::readonly),"Bending moment"))
+ ((Vector3r,pureCreep,Vector3r(0,0,0),(Attr::readonly),"Pure creep curve, used for comparison in calculation."))
+ ((Real,kDam,0,(Attr::readonly),"Damage coefficient on bending, computed from maximum bending moment reached and pure creep behaviour. Its values will vary between :yref:`InelastCohFrictPhys::kr` and :yref:`InelastCohFrictPhys::kRCrp` ."))
// internal attributes
,
createIndex();
);
/// Indexable
- REGISTER_CLASS_INDEX(InelastCohFrictPhys,CohFrictPhys);
+ REGISTER_CLASS_INDEX(InelastCohFrictPhys,FrictPhys);
};
=== modified file 'pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp'
--- pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp 2012-12-14 14:38:44 +0000
+++ pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.cpp 2013-06-25 08:02:13 +0000
@@ -1,5 +1,13 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
#include "Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp"
-
+#include<yade/pkg/dem/ScGeom.hpp>
@@ -11,99 +19,53 @@
InelastCohFrictMat* sdec1 = static_cast<InelastCohFrictMat*>(b1.get());
InelastCohFrictMat* sdec2 = static_cast<InelastCohFrictMat*>(b2.get());
ScGeom6D* geom = YADE_CAST<ScGeom6D*>(interaction->geom.get());
-
- //Create cohesive interractions only once
- if (setCohesionNow && cohesionDefinitionIteration==-1) cohesionDefinitionIteration=scene->iter;
- if (setCohesionNow && cohesionDefinitionIteration!=-1 && cohesionDefinitionIteration!=scene->iter) {
- cohesionDefinitionIteration = -1;
- setCohesionNow = 0;}
+
+ //FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
if (geom) {
if (!interaction->phys) {
interaction->phys = shared_ptr<InelastCohFrictPhys>(new InelastCohFrictPhys());
InelastCohFrictPhys* contactPhysics = YADE_CAST<InelastCohFrictPhys*>(interaction->phys.get());
-// Real Ea = sdec1->young;
-// Real Eb = sdec2->young;
-// Real Va = sdec1->poisson;
-// Real Vb = sdec2->poisson;
- Real Da = geom->radius1;
- Real Db = geom->radius2;
- Real fa = sdec1->frictionAngle;
- Real fb = sdec2->frictionAngle;
- //Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
-
+ Real pi = 3.14159265;
+ Real r1 = geom->radius1;
+ Real r2 = geom->radius2;
+ Real f1 = sdec1->frictionAngle;
+ Real f2 = sdec2->frictionAngle;
+
+ contactPhysics->tangensOfFrictionAngle = tan(min(f1,f2));
+
// harmonic average of modulus
- Real ETC = 2.0*sdec1->eTC*sdec2->eTC/(sdec1->eTC+sdec2->eTC);
- Real ETT = 2.0*sdec1->eTT*sdec2->eTT/(sdec1->eTT+sdec2->eTT);
-
- Real EB = 2.0*sdec1->eB*sdec2->eB/(sdec1->eB+sdec2->eB);
- Real GTw = 2.0*sdec1->gTw*sdec2->gTw/(sdec1->gTw+sdec2->gTw);
-
- contactPhysics->crpT = std::min(sdec1->creepT,sdec2->creepT);
- contactPhysics->crpB = std::min(sdec1->creepB,sdec2->creepB);
- contactPhysics->crpTw = std::min(sdec1->creepTw,sdec2->creepTw);
- //
-
-// Real Ks;
-// if (Va && Vb) Ks = 2.0*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
-// else Ks=0;
-
- Vector3r segment = (Body::byId(interaction->getId1(),scene)->state->pos) - (Body::byId(interaction->getId2(),scene)->state->pos);
- Real length = segment.norm();
- Real pi = 3.14159265;
- Real area = (pow(std::min(Db, Da),2))*pi;
- Real iG = (pow(std::min(2*Db, 2*Da),4))*pi/32.;
- Real iB = (pow(std::min(2*Db, 2*Da),4))*pi/64.;
- // Ignacio Olmedo-Manich, non size-dependent mechanical params.
-
- contactPhysics->knT = ETT*area/length;
- contactPhysics->knC = ETC*area/length;
- contactPhysics->kt = GTw*iG/length;
- contactPhysics->ks = 12*EB*iB/(pow(length,3));
- contactPhysics->kr = EB*iB/length;
-
- contactPhysics->crpT = std::min(sdec1->creepT,sdec2->creepT);
- contactPhysics->crpB = std::min(sdec1->creepB,sdec2->creepB);
- contactPhysics->crpTw = std::min(sdec1->creepTw,sdec2->creepTw);
-
- contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
-
- contactPhysics->maxElastB = iB*std::min(sdec1->sigmaB,sdec2->sigmaB);
- contactPhysics->maxElastTw = iB*std::min(sdec1->sigmaTw,sdec2->sigmaTw);
-
- contactPhysics->unldB = std::min(sdec1->unloadB,sdec2->unloadB);
- contactPhysics->unldT = std::min(sdec1->unloadT,sdec2->unloadT);
- contactPhysics->unldTw = std::min(sdec1->unloadTw,sdec2->unloadTw);
-
- contactPhysics->dElT = std::min(sdec1->disElT,sdec2->disElT);
- contactPhysics->dElC = std::min(sdec1->disElC,sdec2->disElC);
-
- contactPhysics->epsMaxT = std::min(sdec1->epsilonMaxT,sdec2->epsilonMaxT);
- contactPhysics->epsMaxC = std::min(sdec1->epsilonMaxC,sdec2->epsilonMaxC);
-
- contactPhysics->phBMax = std::min(sdec1->phiBMax,sdec2->phiBMax);
- contactPhysics->phTwMax = std::min(sdec1->phiTwMax,sdec2->phiTwMax);
-
- if ((setCohesionOnNewContacts || setCohesionNow) && sdec1->isCohesive && sdec2->isCohesive)
- {
- contactPhysics->cohesionBroken = false;
- contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(Db, Da),2);
- contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(Db, Da),2);
- geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
- }
-
- //contactPhysics->elasticRollingLimit = elasticRollingLimit;
- }
- else {// !isNew, but if setCohesionNow, all contacts are initialized like if they were newly created
- InelastCohFrictPhys* contactPhysics = YADE_CAST<InelastCohFrictPhys*>(interaction->phys.get());
- if ((setCohesionNow && sdec1->isCohesive && sdec2->isCohesive) || contactPhysics->initCohesion)
- {
- contactPhysics->cohesionBroken = false;
- contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
- contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
- geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
- contactPhysics->initCohesion=false;
- }
+ contactPhysics->knC = 2.0*sdec1->compressionModulus*r1*sdec2->compressionModulus*r2/(sdec1->compressionModulus*r1+sdec2->compressionModulus*r2);
+ contactPhysics->knT = 2.0*sdec1->tensionModulus*r1*sdec2->tensionModulus*r2/(sdec1->tensionModulus*r1+sdec2->tensionModulus*r2);
+ contactPhysics->ks = 2.0*sdec1->shearModulus*r1*sdec2->shearModulus*r2/(sdec1->shearModulus*r1+sdec2->shearModulus*r2);
+
+ // harmonic average of coeficients for bending and twist coeficients
+ Real AlphaKr = 2.0*sdec1->alphaKr*sdec2->alphaKr/(sdec1->alphaKr+sdec2->alphaKr);
+ Real AlphaKtw = 2.0*sdec1->alphaKtw*sdec2->alphaKtw/(sdec1->alphaKtw+sdec2->alphaKtw);
+
+ contactPhysics->kr = r1*r2*contactPhysics->ks*AlphaKr;
+ contactPhysics->ktw = r1*r2*contactPhysics->ks*AlphaKtw;
+
+ contactPhysics->kTCrp = contactPhysics->knT*min(sdec1->creepTension,sdec2->creepTension);
+ contactPhysics->kRCrp = contactPhysics->kr*min(sdec1->creepBending,sdec2->creepBending);
+ contactPhysics->kTwCrp = contactPhysics->ktw*min(sdec1->creepTwist,sdec2->creepTwist);
+
+ contactPhysics->kRUnld = contactPhysics->kr*min(sdec1->unloadBending,sdec2->unloadBending);
+ contactPhysics->kTUnld = contactPhysics->knT*min(sdec1->unloadTension,sdec2->unloadTension);
+ contactPhysics->kTwUnld = contactPhysics->ktw*min(sdec1->unloadTwist,sdec2->unloadTwist);
+
+ contactPhysics->maxElC = min(sdec1->sigmaCompression,sdec2->sigmaCompression)*pow(min(r2, r1),2);
+ contactPhysics->maxElT = min(sdec1->sigmaTension,sdec2->sigmaTension)*pow(min(r2, r1),2);
+ contactPhysics->maxElB = min(sdec1->nuBending,sdec2->nuBending)*pow(min(r2, r1),3);
+ contactPhysics->maxElTw = min(sdec1->nuTwist,sdec2->nuTwist)*pow(min(r2, r1),3);
+
+ contactPhysics->shearAdhesion = min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(min(r1, r2),2);
+
+ contactPhysics->maxExten = min(sdec1->epsilonMaxTension*r1,sdec2->epsilonMaxTension*r2);
+ contactPhysics->maxContract = min(sdec1->epsilonMaxCompression*r1,sdec2->epsilonMaxCompression*r2);
+
+ contactPhysics->maxBendMom = min(sdec1->etaMaxBending,sdec2->etaMaxBending)*pow(min(r2, r1),3);
+ contactPhysics->maxTwist = 2*pi*min(sdec1->etaMaxTwist,sdec2->etaMaxTwist);
}
}
};
=== modified file 'pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp'
--- pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp 2012-12-14 14:38:44 +0000
+++ pkg/dem/Ip2_2xInelastCohFrictMat_InelastCohFrictPhys.hpp 2013-06-25 08:02:13 +0000
@@ -1,3 +1,11 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
#pragma once
#include "Ip2_CohFrictMat_CohFrictMat_CohFrictPhys.hpp"
#include "InelastCohFrictMat.hpp"
@@ -15,8 +23,6 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xInelastCohFrictMat_InelastCohFrictPhys,IPhysFunctor,
"Generates cohesive-frictional interactions with moments. Used in the contact law :yref:`Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment`.",
- ((bool,setCohesionNow,false,,"If true, assign cohesion to all existing contacts in current time-step. The flag is turned false automatically, so that assignment is done in the current timestep only."))
- ((bool,setCohesionOnNewContacts,false,,"If true, assign cohesion at all new contacts. If false, only existing contacts can be cohesive (also see :yref:`Ip2_CohFrictMat_CohFrictMat_CohFrictPhys::setCohesionNow`), and new contacts are only frictional."))
,
cohesionDefinitionIteration = -1;
);
=== modified file 'pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp'
--- pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp 2013-05-14 21:24:11 +0000
+++ pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.cpp 2013-06-25 08:02:13 +0000
@@ -1,7 +1,17 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
+
+
#include "Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp"
Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy()
-{
+{ //FIXME : this have to be checked and adapted
Real normEnergy=0;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
@@ -13,7 +23,7 @@
return normEnergy;
}
Real Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy()
-{
+{ //FIXME : this have to be checked and adapted
Real shearEnergy=0;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
@@ -29,251 +39,188 @@
void Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
{
-// cout<<"Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment"<<endl;
- const Real& dt = scene->dt;
+//FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
const int &id1 = contact->getId1();
const int &id2 = contact->getId2();
+ const Real& dt = scene->dt;
ScGeom6D* geom = YADE_CAST<ScGeom6D*> (ig.get());
InelastCohFrictPhys* phys = YADE_CAST<InelastCohFrictPhys*> (ip.get());
- Vector3r& shearForce = phys->shearForce;
-
- cout<<"id1= "<<id1<<" id2= "<<id2<<endl;
-
- if (contact->isFresh(scene)) shearForce = Vector3r::Zero();
-
- ///Tension-Compresion///
-
- Real un = geom->penetrationDepth;
- Real Fn = phys->knT*(un);
- //FIXME: Check signs on TESTS
+ if (contact->isFresh(scene)) phys->shearForce = Vector3r::Zero();
+
+ Real un = geom->penetrationDepth-phys->unp;
+ Real Fn;
+
State* de1 = Body::byId(id1,scene)->state.get();
State* de2 = Body::byId(id2,scene)->state.get();
- if(!(phys->isBrokenT)){
-
- /// Tension ///
- if(un<=0){
- if((un>= -(phys->dElT)) && !onplastT){
- Fn = phys->knT*(un);
-// cout<<"TENSION ELASTIC Fn= "<<Fn<<" un= "<<un<<endl;
- }
- if(un< -(phys->dElT) || onplastT){
- onplastT = true;
- Fn = phys->knT*(-phys->dElT) + phys->crpT*(un+phys->dElT);
-// cout<<"TENSION PLASTIC Fn= "<<Fn<<" un= "<<un<<" normalF= "<<phys->normalForce.norm()<<endl;
-
- if(phys->unloadedT ||(-phys->normalForce.norm()<=Fn)){
- Fn = phys->knT*(-phys->dElT)+ phys->crpT*(lastPlastUn+phys->dElT) + phys->unldT*(un-lastPlastUn);
- phys->unloadedT = true;
-// cout<<"TENSION PLASTIC unloading Fn= "<<Fn<<" un= "<<un<<endl;
-
- if(un<=lastPlastUn){ // Recovers creep after unload and reload
- phys->unloadedT=false;
- Fn = phys->knT*(-phys->dElT) + phys->crpT*(un+phys->dElT);
-// cout<<"TENSION PLASTIC load aftr unld Fn= "<<Fn<<" un= "<<un<<" LPun= "<<lastPlastUn<<endl;
- }
-
- } // Unloading:: //FIXME: ?? Check Fn Sign
- else{// loading, applying Creeping
- lastPlastUn = un;
-// cout<<"TENSION PLASTIC Creep Fn= "<<Fn<<" un= "<<un<<endl;
-
- if (un<-phys->epsMaxT){ // Plastic rupture //
- Fn = 0.0;
- phys->isBrokenT = true;
-// cout<<"TENSION PLASTIC creep BROKEN Fn= "<<Fn<<" un= "<<un<<endl;
- }
- }
- }
- }
- /// Compresion ///
- if(un>0){
- if((un<=phys->dElC) && !onplastT){
- Fn = phys->knC*(un);
-// cout<<"COMPRESION ELASTIC Fn= "<<Fn<<" un= "<<un<<endl;
-
- }
- if(un>phys->dElC || onplastT){
- onplastT = true;
- Fn = phys->knC*(phys->dElC) + phys->crpT*(un-phys->dElC);
- if(phys->unloadedC || (phys->normalForce.norm()>=Fn)){
- Fn = phys->knC*(phys->dElC)+ phys->crpT*(lastPlastUn-phys->dElC) + phys->unldT*(un-lastPlastUn);
- phys->unloadedC = true;
-// cout<<"COMPRESION PLASTIC unloading Fn= "<<Fn<<" un= "<<un<<endl;
-
- if(un>=lastPlastUn){ // Recovers creep after unload and reload
- phys->unloadedC=false;
- Fn = phys->knC*(phys->dElC) + phys->crpT*(un-phys->dElC);
-// cout<<"COMPRESION PLASTIC load aftr unld Fn= "<<Fn<<" un= "<<un<<" LPun= "<<lastPlastUn<<endl;
- }
-
- } // Unloading:: //FIXME: Verify Fn Sign
- else{// loading, applying Creeping
- lastPlastUn = un;
-// cout<<"COMPRESION PLASTIC Creep Fn= "<<Fn<<" un= "<<un<<endl;
- // Fn stills Fn = phys->crpT*(un-phys->unp) ;}
- if (un>phys->epsMaxC){ // Plastic rupture //
- Fn = 0.0;
- phys->isBrokenT = true;
-// cout<<"COMPRESION PLASTIC creep BROKEN Fn= "<<Fn<<" un= "<<un<<endl;
- }
- }
- }
- }
- phys->normalForce = Fn*geom->normal;
- }
-// if ((-Fn)> phys->normalAdhesion) {//normal plasticity
-// Fn=-phys->normalAdhesion;
-// phys->unp = un+phys->normalAdhesion/phys->kn;
-// if (phys->unpMax && phys->unp<phys->unpMax)
-// scene->interactions->requestErase(contact); return;
-// }
-
- // cout<<"Tension-Comp normalF= "<<phys->normalForce<<endl;
- ///end tension-compression///
-
- ///Shear/// ELASTOPLASTIC perfect law TO BE DONE
- //FIXME:: TO DO::Shear ElastoPlastic perfect LAW
- ///////////////////////// CREEP START ///////////
- if (shear_creep) shearForce -= phys->ks*(shearForce*dt/creep_viscosity);
- ///////////////////////// CREEP END ////////////
-
- shearForce = geom->rotate(phys->shearForce);
+ if(un<=0){/// tension ///
+ if(-un>phys->maxExten || phys->isBroken){//plastic failure.
+ phys->isBroken=1;
+ phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
+ scene->interactions->requestErase(contact);
+ return;
+ }
+ Fn=phys->knT*un; //elasticity
+ if(-Fn>phys->maxElT || phys->onPlastT){ //so we are on plastic deformation.
+ phys->onPlastT=1;
+ phys->onPlastC=1; //if plasticity is reached on tension, set it to compression too.
+ if(phys->maxCrpRchdT[0]<un){ //unloading/reloading on plastic deformation.
+ Fn = phys->kTUnld*(un-phys->maxCrpRchdT[0])+phys->maxCrpRchdT[1];
+ }
+ else{//loading on plastic deformation : creep.
+ Fn = -phys->maxElT+phys->kTCrp*(un+phys->maxElT/phys->knT);
+ phys->maxCrpRchdT[0]=un; //new maximum is reached.
+ phys->maxCrpRchdT[1]=Fn;
+ }
+ if (Fn>0){ //so the contact just passed the equilibrium state, set new "unp" who stores the plastic equilibrium state.
+ phys->unp=geom->penetrationDepth;
+ phys->maxCrpRchdT[0]=1e20;
+ phys->maxElT=0;
+ }
+ }
+ else{ //elasticity
+ phys->maxCrpRchdT[0]=un;
+ phys->maxCrpRchdT[1]=Fn;
+ }
+ }
+
+ else{/// compression /// similar to tension.
+ if(un>phys->maxContract || phys->isBroken){
+ phys->isBroken=1;
+ phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
+ if(geom->penetrationDepth<=0){ //do not erase the contact while penetrationDepth<0 because it would be recreated at next timestep.
+ scene->interactions->requestErase(contact);
+ }
+ return;
+ }
+ Fn=phys->knC*un;
+ if(Fn>phys->maxElC || phys->onPlastC){
+ phys->onPlastC=1;
+ if(phys->maxCrpRchdC[0]>un){
+ Fn = phys->kTUnld*(un-phys->maxCrpRchdC[0])+phys->maxCrpRchdC[1];
+ }
+ else{
+ Fn = phys->maxElC+phys->kTCrp*(un-phys->maxElC/phys->knC);
+ phys->maxCrpRchdC[0]=un;
+ phys->maxCrpRchdC[1]=Fn;
+ }
+ if (Fn<0){
+ phys->unp=geom->penetrationDepth;
+ phys->maxCrpRchdC[0]=-1e20;
+ phys->maxElC=0;
+ }
+ }
+ else{
+ phys->maxCrpRchdC[0]=un;
+ phys->maxCrpRchdC[1]=Fn;
+ }
+ }
+
+ /// Shear ///
+ Vector3r shearForce = geom->rotate(phys->shearForce);
const Vector3r& dus = geom->shearIncrement();
//Linear elasticity giving "trial" shear force
- shearForce -= phys->ks*dus;
-
-
- Real Fs = phys->shearForce.norm();
+ shearForce += phys->ks*dus;
+ Real Fs = shearForce.norm();
Real maxFs = phys->shearAdhesion;
-
-
- if (!phys->cohesionDisablesFriction || maxFs==0)
- maxFs += Fn*phys->tangensOfFrictionAngle;
- maxFs = std::max((Real) 0, maxFs);
+ if (maxFs==0)maxFs = Fn*phys->tangensOfFrictionAngle;
+ maxFs = std::max((Real) 0, maxFs);
if (Fs > maxFs) {//Plasticity condition on shear force
-// cout<<"Plastshear ShearAdh= "<<phys->shearAdhesion<<endl;
- if (phys->fragile && !phys->cohesionBroken) {
- phys->SetBreakingState();
+ if (!phys->cohesionBroken) {
+ phys->cohesionBroken=1;
+ phys->shearAdhesion=0;
maxFs = max((Real) 0, Fn*phys->tangensOfFrictionAngle);
}
maxFs = maxFs / Fs;
- Vector3r trialForce=shearForce;
shearForce *= maxFs;
- if (scene->trackEnergy){
- Real dissip=((1/phys->ks)*(trialForce-shearForce))/*plastic disp*/ .dot(shearForce)/*active force*/;
- if(dissip>0) scene->energy->add(dissip,"plastDissip",plastDissipIx,/*reset*/false);}
- if (Fn<0) phys->normalForce = Vector3r::Zero();//Vector3r::Zero()
- }
-// cout<<"Fs= "<<Fs<<" maxFs= "<<maxFs<<endl;
-
- //Apply the force
- applyForceAtContactPoint(-phys->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
-// Vector3r force = -phys->normalForce-shearForce;
-// scene->forces.addForce(id1,force);
-// scene->forces.addForce(id2,-force);
-// scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
-// scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
- /// end Shear ///
-
- /// Moment law ///
- /// Bending///
- if(!(phys->isBrokenB)){
- Vector3r relAngVel = geom->getRelAngVel(de1,de2,dt);
- Vector3r relAngVelBend = relAngVel - geom->normal.dot(relAngVel)*geom->normal; // keep only the bending part
- Vector3r relRotBend = relAngVelBend*dt; // relative rotation due to rolling behaviour
- Vector3r& momentBend = phys->moment_bending;
- momentBend = geom->rotate(momentBend); // rotate moment vector (updated)
-
- //To check if in elastic zone in current iteration
- Real BendValue = (phys->moment_bending-phys->kr*relRotBend).norm();
- Real MaxElastB = phys->maxElastB;
- bool elasticBState = (BendValue<=MaxElastB);
-
- if(!onplastB && elasticBState){
- momentBend = momentBend-phys->kr*relRotBend;
-// cout<<"BENDING Elastic"<<" momentB= "<<momentBend<<endl;
-
- }else{ ///Bending Plasticity///
- onplastB = true;
- BendValue = (phys->moment_bending-phys->crpB*relRotBend).norm();
-// cout<<"BENDING Plastic"<<" momentB= "<<momentBend<<endl;
- // Unloading:: RelRotBend > 0 ::::
- if(phys->unloadedB || phys->moment_bending.norm()>=BendValue){
- momentBend = momentBend-phys->unldB*relRotBend;
- phys->unloadedB = true;
-// cout<<"BENDING Plastic UNLD"<<" momentB= "<<momentBend<<endl;
- if(BendValue>=lastPlastBend){phys->unloadedB = false;}
-
- }
- else{
- momentBend = momentBend-phys->crpB*relRotBend;// loading, applying Creeping
- Vector3r AbsRot = de1->rot()-de2->rot();
- Real AbsBending = (AbsRot - geom->normal.dot(AbsRot)*geom->normal).norm();
- lastPlastBend= BendValue;
- if (AbsBending>phys->phBMax){ // Plastic rupture //FIXME: This line is not ok, need to find the compare the total(or just plastic) angular displacement to PhBmax, true meaning of relRotBend??
- momentBend = Vector3r(0,0,0);
- phys->isBrokenB = true;
-// cout<<"BENDING Plastic BREAK"<<" momentB= "<<momentBend<<endl;
- }
- }
- }
- phys->moment_bending = momentBend;
- }
-
- ///Twist///
-
- if(!(phys->isBrokenTw)){
- Vector3r relAngVel = geom->getRelAngVel(de1,de2,dt);
- Vector3r relAngVelTwist = geom->normal.dot(relAngVel)*geom->normal;
- Vector3r relRotTwist = relAngVelTwist*dt; // component of relative rotation along n FIXME: sign?
- Vector3r& momentTwist = phys->moment_twist;
- momentTwist = geom->rotate(momentTwist); // rotate moment vector (updated)
-
- //To check if in elastic zone in current iteration
- Real TwistValue = (phys->moment_twist-phys->kt*relRotTwist).norm();
- Real MaxElastTw = phys->maxElastTw;
- bool elasticTwState = (TwistValue<=MaxElastTw);
-
- if (!onplastTw && elasticTwState){
- momentTwist = momentTwist-phys->kt*relRotTwist; // FIXME: sign?
-// cout<<"TWIST Elast"<<" momentTwist="<<momentTwist<<endl;
- }else { ///Twist Plasticity///
- onplastTw = true;
- TwistValue = (phys->moment_twist-phys->crpTw*relRotTwist).norm();
-// cout<<"TWIST Plast"<<endl;
- if (phys->unloadedTw || phys->moment_twist.norm()>=TwistValue){// Unloading:: RelRotTwist > 0
- momentTwist = momentTwist-phys->unldTw*relRotTwist;
- phys->unloadedTw = true;
-// cout<<"TWIST Plast UNLD"<<endl;
- if(TwistValue>=lastPlastTw){phys->unloadedTw = false;}
- }
-
- else{momentTwist = momentTwist-phys->crpTw*relRotTwist;// loading, applying Creeping
- Vector3r AbsRot = de1->rot()-de2->rot();
- Real AbsTwist = (geom->normal.dot(AbsRot)*geom->normal).norm();
- lastPlastTw= TwistValue;
-// cout<<"TWIST Creep momentTwist="<<momentTwist<<" AbsTwist="<<AbsTwist<<endl;
- if (AbsTwist>phys->phTwMax){ // Plastic rupture //FIXME: This line is not ok, need to find the compare the total(or just plastic) angular displacement to PhBmax, true meaning of relRotBend??
- momentTwist = Vector3r(0,0,0);
- phys->isBrokenTw = true;
- }
- }
- }
- phys->moment_twist = momentTwist;
- }
-
-// cout<<"moment Twist= "<<phys->moment_twist<<endl;
-
-
- // Apply moments now
- Vector3r moment = phys->moment_twist + phys->moment_bending;
- scene->forces.addTorque(id1,-moment);
- scene->forces.addTorque(id2, moment);
- /// Moment law END ///
-
+ }
+
+ //rotational moment are only applied if the cohesion is not broken.
+ /// Twist /// the twist law is driven by twist displacement ("getTwist()").
+ if(!phys->cohesionBroken){
+ Real twist = geom->getTwist() - phys->twp;
+ Real twistM=twist*phys->ktw; //elastic twist moment.
+ bool sgnChanged=0; //whether the twist moment just passed the equilibrium state.
+ if(!contact->isFresh(scene) && phys->moment_twist.dot(twistM*geom->normal)<0)sgnChanged=1;
+ if(abs(twist)>phys->maxTwist){
+ phys->cohesionBroken=1;
+ twistM=0;
+ }
+ else{
+ if(abs(twistM)>phys->maxElTw || phys->onPlastTw){ //plastic deformation.
+ phys->onPlastTw=1;
+ if(abs(phys->maxCrpRchdTw[0])>abs(twist)){ //unloading/reloading
+ twistM = phys->kTwUnld*(twist-phys->maxCrpRchdTw[0])+phys->maxCrpRchdTw[1];
+ }
+ else{//creep loading.
+ int sign = twist<0?-1:1;
+ twistM = sign*phys->maxElTw+phys->kTwCrp*(twist-sign*phys->maxElTw/phys->ktw); //creep
+ phys->maxCrpRchdTw[0]=twist; //new maximum reached
+ phys->maxCrpRchdTw[1]=twistM;
+ }
+ if(sgnChanged){
+ phys->maxElTw=0;
+ phys->twp=geom->getTwist();
+ phys->maxCrpRchdTw[0]=0;
+ }
+ }
+ else{ //elasticity
+ phys->maxCrpRchdTw[0]=twist;
+ phys->maxCrpRchdTw[1]=twistM;
+ }
+ }
+ phys->moment_twist = twistM * geom->normal;
+ }
+ else phys->moment_twist=Vector3r(0,0,0);
+
+ /// Bending /// incremental form.
+ if(!phys->cohesionBroken){
+ Vector3r bendM = phys->moment_bending;
+ Vector3r relAngVel = geom->getRelAngVel(de1,de2,dt);
+ Vector3r relRotBend = (relAngVel - geom->normal.dot(relAngVel)*geom->normal)*dt; // relative rotation due to rolling behaviour
+ bendM = geom->rotate(phys->moment_bending); // rotate moment vector (updated)
+ phys->pureCreep=geom->rotate(phys->pureCreep); // pure creep is updated to compute the damage.
+ Vector3r bendM_elast = bendM-phys->kr*relRotBend;
+ if(bendM_elast.norm()>phys->maxElB || phys->onPlastB){ // plastic behavior
+ phys->onPlastB=1;
+ bendM=bendM-phys->kDam*relRotBend; //trial bending
+ if(bendM.norm()<phys->moment_bending.norm()){ // if bending decreased, we are unloading ...
+ bendM = bendM+phys->kDam*relRotBend-phys->kRUnld*relRotBend; // ... so undo bendM and apply unload coefficient.
+ Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend; // trial pure creep.
+ phys->pureCreep = newPureCreep.norm()<phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while unloading, pure creep must decrease.
+ phys->kDam=phys->kr+(phys->kRCrp-phys->kr)*(phys->maxCrpRchdB.norm()-phys->maxElB)/(phys->maxBendMom-phys->maxElB); // compute the damage coefficient.
+ }
+ else{ // bending increased, so we are loading (bendM has to be unchanged).
+ Vector3r newPureCreep = phys->pureCreep-phys->kRCrp*relRotBend;
+ phys->pureCreep = newPureCreep.norm()>phys->pureCreep.norm()?newPureCreep:phys->pureCreep+phys->kRCrp*relRotBend; // while loading, pure creep must increase.
+ if(phys->pureCreep.norm()<bendM.norm()) bendM=phys->pureCreep; // bending moment can't be greather than pure creep.
+ if(phys->pureCreep.norm()>phys->maxCrpRchdB.norm()) phys->maxCrpRchdB=phys->pureCreep; // maxCrpRchdB must follow the maximum of pure creep.
+ if(phys->pureCreep.norm()>phys->maxBendMom){
+ phys->cohesionBroken=1;
+ bendM=bendM_elast=Vector3r(0,0,0);
+ }
+ }
+ phys->moment_bending=bendM;
+ }
+ else{//elasticity
+ phys->pureCreep=phys->moment_bending=phys->maxCrpRchdB=bendM_elast;
+ phys->kDam=phys->kRCrp;
+ }
+ }
+ phys->shearForce=shearForce;
+ phys->normalForce=-Fn*geom->normal;
+ applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
+ scene->forces.addTorque(id1,-phys->moment_bending-phys->moment_twist);
+ scene->forces.addTorque(id2,phys->moment_bending+phys->moment_twist);
}
YADE_PLUGIN((Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment));
+
+
+
+
+
+
+
=== modified file 'pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp'
--- pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 2013-05-29 09:48:51 +0000
+++ pkg/dem/Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment.hpp 2013-06-25 08:02:13 +0000
@@ -1,3 +1,11 @@
+/*************************************************************************
+* Copyright (C) 2012 by Ignacio Olmedo nolmedo.manich@xxxxxxxxx *
+* Copyright (C) 2012 by François Kneib francois.kneib@xxxxxxxxx *
+* This program is free software; it is licensed under the terms of the *
+* GNU General Public License v2 or later. See file LICENSE for details. *
+*************************************************************************/
+
+
#pragma once
#include "CohesiveFrictionalContactLaw.hpp"
#include "InelastCohFrictPhys.hpp"
@@ -8,20 +16,7 @@
Real shearElastEnergy();
virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. The maximum value of moments can be defined and takes the form of rolling friction. Cohesive -type moment may also be included in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
- ((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
- ((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
- ((bool,shear_creep,false,,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
- ((bool,twist_creep,false,,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
((bool,useIncrementalForm,false,,"use the incremental formulation to compute bending and twisting moments. Creep on the twisting moment is not included in such a case."))
- ((int,plastDissipIx,-1,,"Index for plastic dissipation (with O.trackEnergy)"))
- ((bool,onplastB,false,,"true if plasticity achieved on bending"))
- ((bool,onplastTw,false,,"true if plasticity achieved on twisting"))
- ((bool,onplastT,false,,"true if plasticity achieved on traction/compression")) // if plasticity achieved in traction, also for later compression, and inverse
- ((Real,lastPlastUn,0,,"last tension compression un value of plasticity"))
- ((Real,lastPlastBend,0,,"last bend value of plasticity"))
- ((Real,lastPlastTw,0,,"last bend value of plasticity"))
- ((Real,unPrev,0,,"previous timestep penetrationDepth"))
- ((Real,creep_viscosity,1,,"creep viscosity [Pa.s/m]. probably should be moved to Ip2_CohFrictMat_CohFrictMat_CohFrictPhys."))
,,
.def("normElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
.def("shearElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::shearElastEnergy,"Compute shear elastic energy.")
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2013-05-30 20:19:43 +0000
+++ pkg/dem/Shop.cpp 2013-06-20 14:56:35 +0000
@@ -914,7 +914,7 @@
FOREACH(const shared_ptr<Body>& b,*scene->bodies){
if(b->isClump()){
Clump* clumpSt = YADE_CAST<Clump*>(b->shape.get());
- clumpSt->updateProperties(b, 0);
+ clumpSt->updateProperties(b);
}
}
FOREACH(const shared_ptr<Interaction>& ii, *scene->interactions){
=== modified file 'pkg/dem/WirePM.cpp'
--- pkg/dem/WirePM.cpp 2012-01-23 14:43:54 +0000
+++ pkg/dem/WirePM.cpp 2013-06-25 04:00:41 +0000
@@ -10,6 +10,13 @@
#include<yade/core/Scene.hpp>
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/core/Omega.hpp>
+#include "../../lib/base/Math.hpp"
+
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/triangle_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+
+#include<yade/core/Timing.hpp>
YADE_PLUGIN((WireMat)(WireState)(WirePhys)(Ip2_WireMat_WireMat_WirePhys)(Law2_ScGeom_WirePhys_WirePM));
@@ -22,14 +29,50 @@
//BUG: ????? postLoad is called twice,
LOG_TRACE( "WireMat::postLoad - update material parameters" );
- // compute cross-section area
+ // compute cross-section area for single wire
as = pow(diameter*0.5,2)*Mathr::PI;
+ // check for stress strain curve for single wire
if(strainStressValues.empty()) return; // uninitialized object, don't do nothing at all
if(strainStressValues.size() < 2)
- throw invalid_argument("WireMat.strainStressValues: at least two points must be given.");
+ throw std::invalid_argument("WireMat.strainStressValues: at least two points must be given.");
if(strainStressValues[0](0) == 0. && strainStressValues[0](1) == 0.)
- throw invalid_argument("WireMat.strainStressValues: Definition must start with values greather then zero (strain>0,stress>0)");
+ throw std::invalid_argument("WireMat.strainStressValues: Definition must start with values greater than zero (strain>0,stress>0)");
+
+ switch(type) {
+ case 0:
+ LOG_DEBUG("WireMat - Bertrand's approach");
+ if(!strainStressValuesDT.empty())
+ throw std::invalid_argument("Use of WireMat.strainStressValuesDT has no effect!");
+ break;
+ case 1:
+ // check stress strain curve four double twist if type=1
+ LOG_DEBUG("WireMat - New approach with two curves");
+ if(isDoubleTwist) {
+ if(strainStressValuesDT.empty())
+ throw runtime_error("WireMat.strainStressValuesDT not defined");
+ if(strainStressValuesDT.size() < 2)
+ throw std::invalid_argument("WireMat.strainStressValuesDT: at least two points must be given.");
+ if(strainStressValuesDT[0](0) == 0. && strainStressValuesDT[0](1))
+ throw std::invalid_argument("WireMat.strainStressValuesDT: Definition must start with values greater than zero (strain>0,stress>0)");
+ }
+ break;
+ case 2:
+ // check stress strain curve four double twist if type=2
+ LOG_DEBUG("WireMat - New approach with two curves and initial shift");
+ if(isDoubleTwist) {
+ if(strainStressValuesDT.empty())
+ throw runtime_error("WireMat.strainStressValuesDT not defined");
+ if(strainStressValuesDT.size() < 2)
+ throw std::invalid_argument("WireMat.strainStressValuesDT: at least two points must be given.");
+ if(strainStressValuesDT[0](0) == 0. && strainStressValuesDT[0](1))
+ throw std::invalid_argument("WireMat.strainStressValuesDT: Definition must start with values greater than zero (strain>0,stress>0)");
+ }
+ break;
+ default:
+ throw std::invalid_argument("WireMat.type: Type must be 0, 1 or 2.");
+ break;
+ }
}
@@ -53,6 +96,7 @@
/* get reference to values since values are updated/changed in order to take unloading into account */
vector<Vector2r> &DFValues = phys->displForceValues;
vector<Real> &kValues = phys->stiffnessValues;
+ Real kn = phys->kn;
Real D = displN - phys->initD; // interparticular distance is computed depending on the equilibrium distance
@@ -74,29 +118,29 @@
/* compute normal force Fn */
Real Fn = 0.;
- if ( D > DFValues[0](0) )
- Fn = kValues[0] * (D-phys->plastD);
- else {
- bool isDone = false;
- unsigned int i = 0;
- while (!isDone && i < DFValues.size()) {
- i++;
+ if ( D > DFValues[0](0) ) { // unloading
+ LOG_TRACE("WirePM: Unloading");
+ Fn = kn * (D-phys->plastD);
+ }
+ else { // loading
+ LOG_TRACE("WirePM: Loading");
+ for (unsigned int i=1; i<DFValues.size(); i++) {
if ( D > DFValues[i](0) ) {
- Fn = DFValues[i-1](1) + (D-DFValues[i-1](0))*kValues[i];
- phys->plastD = D - Fn/kValues[0];
- isDone = true;
+ Fn = DFValues[i-1](1) + (D-DFValues[i-1](0))*kValues[i-1];
+ phys->plastD = D - Fn/kn;
// update values for unloading
DFValues[0](0) = D;
DFValues[0](1) = Fn;
+ break;
}
}
}
- TRVAR3( displN, D, Fn );
-
/* compression forces cannot be applied to wires */
if (Fn > 0.) Fn = 0.;
+ TRVAR3( displN, D, Fn );
+
phys->normalForce = Fn*geom->normal; // NOTE: normal is position2-position1 - It is directed from particle1 to particle2
/* compute a limit value to check how far the interaction is from failing */
@@ -118,7 +162,7 @@
}
/* set shear force to zero */
- phys->shearForce = Vector3r::Zero();;
+ phys->shearForce = Vector3r::Zero();
}
@@ -129,6 +173,7 @@
/* avoid any updates if interactions which already exist */
if(interaction->phys) return;
+ //TODO: make boolean to make sure physics are never updated, optimisation of contact detection mesh (no contact detection after link is created)
LOG_TRACE( "Ip2_WireMat_WireMat_WirePhys::go - create interaction physics" );
@@ -137,7 +182,7 @@
/* set equilibrium distance, e.g. initial distance between particle (stress free state) */
shared_ptr<WirePhys> contactPhysics(new WirePhys());
- contactPhysics->initD = geom->penetrationDepth;
+ Real initD = geom->penetrationDepth;
contactPhysics->normalForce = Vector3r::Zero();
/* get values from material */
@@ -151,10 +196,16 @@
if ( mat1->id == mat2->id ) { // interaction of two bodies of the same material
crossSection = mat1->as;
SSValues = mat1->strainStressValues;
- if ( (mat1->isDoubleTwist) && (abs(interaction->getId1()-interaction->getId2())==1) ) // bodies which id differs by 1 are double twisted
+ if ( (mat1->isDoubleTwist) && (abs(interaction->getId1()-interaction->getId2())==1) ) {// bodies which id differs by 1 are double twisted
contactPhysics->isDoubleTwist = true;
- else
+ if ( mat1->type==1 || mat1->type==2 ) {
+ SSValues = mat1->strainStressValuesDT;
+ crossSection *= 2.;
+ }
+ }
+ else {
contactPhysics->isDoubleTwist = false;
+ }
}
else { // interaction of two bodies of two different materials, take weaker material and no double-twist
contactPhysics->isDoubleTwist = false;
@@ -168,28 +219,67 @@
}
}
- if(SSValues.empty()) throw std::invalid_argument("WireMat.strainStressValue is empty!");
-
Real R1 = geom->radius1;
Real R2 = geom->radius2;
- Real l0 = R1 + R2 - contactPhysics->initD; // initial length of the wire (can be single or double twisted)
+ Real l0 = R1 + R2 - initD; // initial length of the wire (can be single or double twisted)
+
+ /* compute displacement-force values */
+ vector<Vector2r> DFValues;
+ vector<Real> kValues;
+ Real dl = 0.;
+ bool isShifted = false;
+
+ /* account for random distortion if type=2 */
+ if ( mat1->type==2 ) {
+ isShifted = true;
+ if (mat1->seed==-1)
+ dl = l0*mat1->lambdau;
+ else {
+ // initialize random number generator
+ static boost::minstd_rand randGen(mat1->seed!=0?mat1->seed:(int)TimingInfo::getNow(true));
+ static boost::variate_generator<boost::minstd_rand&, boost::triangle_distribution<Real> > rnd(randGen, boost::triangle_distribution<Real>(0,0.5,1));
+ Real rndu = rnd();
+ TRVAR1( rndu );
+ dl = l0*mat1->lambdau*rndu;
+ isShifted = true;
+ }
+ }
+ else if ( mat2->type==2 ) {
+ isShifted = true;
+ if (mat2->seed==-1)
+ dl = l0*mat2->lambdau;
+ else {
+ // initialize random number generator
+ static boost::minstd_rand randGen(mat2->seed!=0?mat2->seed:(int)TimingInfo::getNow(true));
+ static boost::variate_generator<boost::minstd_rand&, boost::triangle_distribution<Real> > rnd(randGen, boost::triangle_distribution<Real>(0,0.5,1));
+ Real rndu = rnd();
+ TRVAR1( rndu );
+ dl = l0*mat2->lambdau*rndu;
+ }
+ }
+ contactPhysics->dL=dl;
+ contactPhysics->isShifted=isShifted;
+
+ // update geometry values
+ l0 += dl;
+ contactPhysics->initD = initD;
/* compute threshold displacement-force values (tension negative since ScGem is used!) */
- vector<Vector2r> DFValues;
for ( vector<Vector2r>::iterator it = SSValues.begin(); it != SSValues.end(); it++ ) {
Vector2r values = Vector2r::Zero();
- values(0) = -(*it)(0)*l0;
+// values(0) = -(*it)(0)*l0;
+ values(0) = -(*it)(0)*l0-dl;
values(1) = -(*it)(1)*crossSection;
DFValues.push_back(values);
}
- /* compute elastic stiffness of single wire */
- vector<Real> kValues;
- Real k = DFValues[0](1) / DFValues[0](0);
+ /* compute elastic stiffness for unloading*/
+ Real k = DFValues[0](1) / (DFValues[0](0)+dl);
- /* update values if the interaction is double twisted */
- if ( contactPhysics->isDoubleTwist ) {
+ /* update values if the interaction is a double twist and type=0 */
+ if ( contactPhysics->isDoubleTwist && mat1->type==0 ) {
+ // type=0 (force displacement values are computed by manipulating the values of the single wire by using the parameters lambdak and lambdaEps)
Real alpha = atan( l0 / (3.*Mathr::PI*mat1->diameter) );
Real kh = k * ( l0*mat1->diameter/crossSection ) / ( 48.*cos(alpha) * ( 41./9.*(1.+mat1->poisson) + 17./4.*pow(tan(alpha),2) ) );
k = 2. * ( mat1->lambdak*kh + (1-mat1->lambdak)*k );
@@ -201,12 +291,34 @@
DFValues[i](1) *= mappingF;
}
}
+ else {
+ // type=1 and type=2 (force displacement values have already been computed by given stress-strain curve)
+ }
- /* store displacement-force values in physics */
- contactPhysics->displForceValues = DFValues;
+ /* store elastic/unloading stiffness as kn in physics */
+ contactPhysics->kn = k;
+ contactPhysics->ks = 0.;
+ TRVAR1( k );
+
+ /* consider an additional point for the initial shift if type==2 */
+ if ( mat1->type==2 ) {
+ Vector2r values = Vector2r::Zero();
+ values(0) = -dl+mat1->lambdaF*(DFValues[0](0)+dl);
+ values(1) = DFValues[0](1)*mat1->lambdaF;
+ k = values(1) / values(0);
+ if ( mat1->lambdaF<1. )
+ DFValues.insert( DFValues.begin(), values );
+ }
+ else if ( mat2->type==2 ) {
+ Vector2r values = Vector2r::Zero();
+ values(0) = -dl+mat2->lambdaF*(DFValues[0](0)+dl);
+ values(1) = DFValues[0](1)*mat2->lambdaF;
+ k = values(1) / values(0);
+ if ( mat2->lambdaF<1. )
+ DFValues.insert( DFValues.begin(), values );
+ }
/* compute stiffness-values of wire */
- contactPhysics->kn = k;
kValues.push_back(k);
for( unsigned int i = 1 ; i < DFValues.size(); i++ ) {
Real deltau = -DFValues[i](0) + DFValues[i-1](0);
@@ -214,6 +326,12 @@
k = deltaF/deltau;
kValues.push_back(k);
}
+
+ /* add zero values for first point */
+ DFValues.insert( DFValues.begin(), Vector2r::Zero() );
+
+ /* store values in physics */
+ contactPhysics->displForceValues = DFValues;
contactPhysics->stiffnessValues = kValues;
/* set particles as linked */
@@ -227,3 +345,4 @@
}
WirePhys::~WirePhys(){}
+
=== modified file 'pkg/dem/WirePM.hpp'
--- pkg/dem/WirePM.hpp 2012-01-15 21:00:38 +0000
+++ pkg/dem/WirePM.hpp 2013-06-25 04:00:41 +0000
@@ -9,7 +9,7 @@
/**
=== OVERVIEW OF WirePM ===
-A particle model to simulate single wires and rockfall meshes (see Bertrad et al. 2005, Bertrad et al. 2008).
+A particle model to simulate single wires and rockfall meshes (see Bertrad et al. 2005, Bertrad et al. 2008, Thoeni et al. 2013).
Features of the interaction law:
@@ -19,8 +19,8 @@
3. The force displacement curve which defines the interaction forces is piecewise linear and defined by the stress-strain curve of the wire material. Any piecewise linear curve can be used.
-Remarks:
-- The contact law is new and still needs some testing :-)
+4. Three different types of wire models are available.
+
*/
#pragma once
@@ -48,12 +48,22 @@
virtual bool stateTypeOk(State* s) const { return (bool)dynamic_cast<WireState*>(s); }
void postLoad(WireMat&);
YADE_CLASS_BASE_DOC_ATTRS_CTOR(WireMat,FrictMat,"Material for use with the Wire classes",
- ((Real,diameter,0.0027,,"Diameter of the single wire in [m] (the diameter is used to compute the cross-section area of the wire)."))
- ((vector<Vector2r>,strainStressValues,,Attr::triggerPostLoad,"Piecewise linear definition of the stress-strain curve by set of points (strain[-]>0,stress[Pa]>0) for one single wire. Tension only is considered and the point (0,0) is not needed!"))
- ((bool,isDoubleTwist,false,,"Type of the mesh. If true two particles of the same material which body ids differ by one will be considered as double-twisted interaction."))
- ((Real,lambdaEps,0.4,,"Parameter between 0 and 1 to reduce the failure strain of the double-twisted wire (as used by [Bertrand2008]_). [-]"))
- ((Real,lambdak,0.21,,"Parameter between 0 and 1 to compute the elastic stiffness of the double-twisted wire (as used by [Bertrand2008]_): $k^D=2(\\lambda_k k_h + (1-\\lambda_k)k^S)$. [-]"))
- ((Real,as,0,Attr::readonly,"Cross-section area of a single wire used for the computation of the limit normal contact forces. [m²]"))
+ ((Real,diameter,0.0027,,"Diameter of the single wire in [m] (the diameter is used to compute the cross-section area of the wire)."))
+ ((unsigned int,type,0,,"Three different types are considered:\n\n"
+ "== ===============================================================\n"
+ "type=0 Corresponds to Bertrand's approach (see [Bertrand2008]_): only one stress-strain curve is used\n"
+ "type=1 New approach: two stress-strain curves (see [Thoeni2013]_)\n"
+ "type=2 New approach: two stress-strain curves with changed initial stiffness and horizontal shift (shift is random if $\\text{seed}\\geq0$, for more details see [Thoeni2013]_)\n"
+ "== ==============================================================="))
+ ((vector<Vector2r>,strainStressValues,,Attr::triggerPostLoad,"Piecewise linear definition of the stress-strain curve by set of points (strain[-]>0,stress[Pa]>0) for one single wire. Tension only is considered and the point (0,0) is not needed! NOTE: Vector needs to be initialized!"))
+ ((vector<Vector2r>,strainStressValuesDT,,Attr::triggerPostLoad,"Piecewise linear definition of the stress-strain curve by set of points (strain[-]>0,stress[Pa]>0) for the double twist. Tension only is considered and the point (0,0) is not needed! If this value is given the calculation will be based on two different stress-strain curves without considering the parameter introduced by [Bertrand2008]_ (see [Thoeni2013]_)."))
+ ((bool,isDoubleTwist,false,,"Type of the mesh. If true two particles of the same material which body ids differ by one will be considered as double-twisted interaction."))
+ ((Real,lambdaEps,0.47,,"Parameter between 0 and 1 to reduce strain at failure of a double-twisted wire (as used by [Bertrand2008]_). [-]"))
+ ((Real,lambdak,0.73,,"Parameter between 0 and 1 to compute the elastic stiffness of a double-twisted wire (as used by [Bertrand2008]_): $k^D=2(\\lambda_k k_h + (1-\\lambda_k)k^S)$. [-]"))
+ ((int,seed,12345,,"Integer used to initialize the random number generator for the calculation of the distortion. If the integer is equal to 0 a internal seed number based on the time is computed. [-]"))
+ ((Real,lambdau,0.2,,"Parameter between 0 and 1 introduced by [Thoeni2013]_ which defines the maximum shift of the force-displacement curve in order to take an additional initial elongation (e.g. wire distortion/imperfections, slipping, system flexibility) into account: $\\Delta l^*=\\lambda_u l_0 \\text{rnd(seed)}$. [-]"))
+ ((Real,lambdaF,1.0,,"Parameter between 0 and 1 introduced by [Thoeni2013]_ which defines where the shifted force-displacement curve intersects with the new initial stiffness: $F^*=\\lambda_F F_{\\text{elastic}}$. [-]"))
+ ((Real,as,0.,Attr::readonly,"Cross-section area of a single wire used to transform stress into force. [m²]"))
,
createIndex();
);
@@ -69,13 +79,15 @@
virtual ~WirePhys();
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(WirePhys,FrictPhys,"Representation of a single interaction of the WirePM type, storage for relevant parameters",
- ((Real,initD,0,,"Equilibrium distance for particles. Computed as the initial inter-particular distance when particle are linked."))
+ ((Real,initD,0.,,"Equilibrium distance for particles. Computed as the initial inter-particular distance when particle are linked."))
((bool,isLinked,false,,"If true particles are linked and will interact. Interactions are linked automatically by the definition of the corresponding interaction radius. The value is false if the wire breaks (no more interaction)."))
((bool,isDoubleTwist,false,,"If true the properties of the interaction will be defined as a double-twisted wire."))
((vector<Vector2r>,displForceValues,,Attr::readonly,"Defines the values for force-displacement curve."))
- ((vector<Real>,stiffnessValues,,Attr::readonly,"Defines the values for the different stiffness (first value corresponds to elastic stiffness kn)."))
- ((Real,plastD,0,Attr::readonly,"Plastic part of the inter-particular distance of the previous step. \n\n.. note::\n\t Only elastic displacements are reversible (the elastic stiffness is used for unloading) and compressive forces are inadmissible. The compressive stiffness is assumed to be equal to zero (see [Bertrand2005]_).\n\n.."))
- ((Real,limitFactor,0.,Attr::readonly,"This value indicates on how far from failing the wire is, e.g. actual normal displacement divided by admissible normal displacement multiplied by actual normal force divided by admissible normal force."))
+ ((vector<Real>,stiffnessValues,,Attr::readonly,"Defines the values for the various stiffnesses (the elastic stiffness is stored as kn)."))
+ ((Real,plastD,0.,Attr::readonly,"Plastic part of the inter-particular distance of the previous step. \n\n.. note::\n\t Only elastic displacements are reversible (the elastic stiffness is used for unloading) and compressive forces are inadmissible. The compressive stiffness is assumed to be equal to zero.\n\n.."))
+ ((Real,limitFactor,0.,Attr::readonly,"This value indicates on how far from failing the wire is, e.g. actual normal displacement divided by admissible normal displacement."))
+ ((bool,isShifted,false,Attr::readonly,"If true :yref:`WireMat` type=2 and the force-displacement curve will be shifted."))
+ ((Real,dL,0.,Attr::readonly,"Additional wire length for considering the distortion for :yref:`WireMat` type=2 (see [Thoeni2013]_)."))
,
createIndex();
,
=== modified file 'py/CMakeLists.txt'
--- py/CMakeLists.txt 2013-03-18 18:47:34 +0000
+++ py/CMakeLists.txt 2013-07-01 18:57:33 +0000
@@ -32,18 +32,18 @@
ADD_LIBRARY(_eudoxos SHARED "${CMAKE_CURRENT_SOURCE_DIR}/_eudoxos.cpp")
SET_TARGET_PROPERTIES(_eudoxos PROPERTIES PREFIX "")
-TARGET_LINK_LIBRARIES(_eudoxos plugins support core)
+TARGET_LINK_LIBRARIES(_eudoxos yade)
INSTALL(TARGETS _eudoxos DESTINATION "${YADE_PY_PATH}/yade/")
ADD_LIBRARY(_utils SHARED "${CMAKE_CURRENT_SOURCE_DIR}/_utils.cpp")
SET_TARGET_PROPERTIES(_utils PROPERTIES PREFIX "")
-TARGET_LINK_LIBRARIES(_utils plugins support core)
+TARGET_LINK_LIBRARIES(_utils yade)
INSTALL(TARGETS _utils DESTINATION "${YADE_PY_PATH}/yade/")
ADD_LIBRARY(_packPredicates SHARED "${CMAKE_CURRENT_SOURCE_DIR}/pack/_packPredicates.cpp")
SET_TARGET_PROPERTIES(_packPredicates PROPERTIES PREFIX "")
-TARGET_LINK_LIBRARIES(_packPredicates plugins support core)
+TARGET_LINK_LIBRARIES(_packPredicates yade)
IF(ENABLE_GTS AND NOT(PY_gts))
TARGET_LINK_LIBRARIES(_packPredicates _gts)
ENDIF(ENABLE_GTS AND NOT(PY_gts))
@@ -53,7 +53,7 @@
ADD_LIBRARY(_packSpheres SHARED "${CMAKE_CURRENT_SOURCE_DIR}/pack/_packSpheres.cpp")
SET_TARGET_PROPERTIES(_packSpheres PROPERTIES PREFIX "")
-TARGET_LINK_LIBRARIES(_packSpheres plugins support core)
+TARGET_LINK_LIBRARIES(_packSpheres yade)
INSTALL(TARGETS _packSpheres DESTINATION "${YADE_PY_PATH}/yade/")
ADD_LIBRARY(_packObb SHARED "${CMAKE_CURRENT_SOURCE_DIR}/pack/_packObb.cpp")
@@ -62,14 +62,10 @@
ADD_LIBRARY(wrapper SHARED "${CMAKE_CURRENT_SOURCE_DIR}/wrapper/yadeWrapper.cpp")
SET_TARGET_PROPERTIES(wrapper PROPERTIES PREFIX "")
-TARGET_LINK_LIBRARIES(wrapper plugins support)
+TARGET_LINK_LIBRARIES(wrapper yade)
INSTALL(TARGETS wrapper DESTINATION "${YADE_PY_PATH}/yade/")
ADD_LIBRARY(_customConverters SHARED "${CMAKE_CURRENT_SOURCE_DIR}/wrapper/customConverters.cpp")
SET_TARGET_PROPERTIES(_customConverters PROPERTIES PREFIX "")
-TARGET_LINK_LIBRARIES(_customConverters plugins support)
+TARGET_LINK_LIBRARIES(_customConverters yade)
INSTALL(TARGETS _customConverters DESTINATION "${YADE_PY_PATH}/yade/")
-
-
-
-
=== modified file 'py/wrapper/yadeWrapper.cpp'
--- py/wrapper/yadeWrapper.cpp 2013-04-17 07:39:24 +0000
+++ py/wrapper/yadeWrapper.cpp 2013-06-21 06:34:26 +0000
@@ -92,7 +92,7 @@
private:
void checkClump(shared_ptr<Body> b){
if (!(b->isClump())){
- PyErr_SetString(PyExc_TypeError,("Error: Body"+lexical_cast<string>(b->getId())+" is not a clump.").c_str());
+ PyErr_SetString(PyExc_TypeError,("Error: Body"+lexical_cast<string>(b->getId())+" is not a clump.").c_str());
python::throw_error_already_set();
}
}
@@ -149,7 +149,7 @@
};
FOREACH(Body::id_t id, ids) Clump::add(clumpBody,Body::byId(id,scene));
- Clump::updateProperties(clumpBody,/*intersecting*/ false);
+ Clump::updateProperties(clumpBody);
return clumpBody->getId();
}
python::tuple appendClump(vector<shared_ptr<Body> > bb){
@@ -166,7 +166,7 @@
if (bp->isClump()){
if (bp == clp) {PyErr_Warn(PyExc_UserWarning,("Warning: Body "+lexical_cast<string>(bid)+" and clump "+lexical_cast<string>(cid)+" are the same bodies. Body was not added.").c_str()); return;}
Clump::add(clp,bp);//add clump bid to clump cid
- Clump::updateProperties(clp,/*intersecting*/false);
+ Clump::updateProperties(clp);
proxee->erase(bid);//erase old clump
return;
}
@@ -175,11 +175,11 @@
shared_ptr<Body> bpClumpPointer = Body::byId(bpClumpId,scene);
if (bpClumpPointer == clp) {PyErr_Warn(PyExc_UserWarning,("Warning: Body "+lexical_cast<string>(bid)+" is already a clump member of clump "+lexical_cast<string>(cid)+". Body was not added.").c_str()); return;}
Clump::add(clp,bpClumpPointer);//add clump bpClumpId to clump cid
- Clump::updateProperties(clp,/*intersecting*/false);
+ Clump::updateProperties(clp);
proxee->erase(bpClumpId);//erase old clump
return;
}
- else {Clump::add(clp,bp); Clump::updateProperties(clp,/*intersecting*/false);}// bp must be a standalone!
+ else {Clump::add(clp,bp); Clump::updateProperties(clp);}// bp must be a standalone!
}
void releaseFromClump(Body::id_t bid, Body::id_t cid){
Scene* scene(Omega::instance().getScene().get()); // get scene
@@ -193,7 +193,7 @@
std::map<Body::id_t,Se3r>& members = clump->members;
if (members.size() == 2) {PyErr_Warn(PyExc_UserWarning,("Warning: Body "+lexical_cast<string>(bid)+" not released from clump "+lexical_cast<string>(cid)+", because number of clump members would get < 2!").c_str()); return;}
Clump::del(clp,bp);//release bid from cid
- Clump::updateProperties(clp,/*intersecting*/false);
+ Clump::updateProperties(clp);
} else { PyErr_Warn(PyExc_UserWarning,("Warning: Body "+lexical_cast<string>(bid)+" must be a clump member of clump "+lexical_cast<string>(cid)+". Body was not released.").c_str()); return;}
} else { PyErr_Warn(PyExc_UserWarning,("Warning: Body "+lexical_cast<string>(bid)+" is not a clump member. Body was not released.").c_str()); return;}
}
@@ -343,105 +343,6 @@
}
return ret;
}
- Real getClumpVolume(const shared_ptr<Body>& b, int numTries){//algorithm for estimation of clump volumes
- checkClump(b);
- Scene* scene(Omega::instance().getScene().get()); // get scene
- shared_ptr<Sphere> sph (new Sphere);
- int Sph_Index = sph->getClassIndexStatic(); // get sphere index for checking if bodies are spheres
- //get boundaries of clump:
- const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(b->shape);
- std::map<Body::id_t,Se3r>& members = clump->members;
- Real xmax = -1e10;
- Real ymax = xmax;
- Real zmax = xmax;
- Real xmin = 1e10;
- Real ymin = xmin;
- Real zmin = xmin;
- FOREACH(MemberMap::value_type& mm, members){
- const Body::id_t& memberId=mm.first;
- const shared_ptr<Body>& member=Body::byId(memberId,scene);
- assert(member->isClumpMember());
- if (member->shape->getClassIndex() == Sph_Index){//clump member should be a sphere
- const Sphere* sphere = YADE_CAST<Sphere*> (member->shape.get());
- xmax = max(member->state->pos[0] + sphere->radius, xmax);
- ymax = max(member->state->pos[1] + sphere->radius, ymax);
- zmax = max(member->state->pos[2] + sphere->radius, zmax);
- xmin = min(member->state->pos[0] - sphere->radius, xmin);
- ymin = min(member->state->pos[1] - sphere->radius, ymin);
- zmin = min(member->state->pos[2] - sphere->radius, zmin);
- }
- else {PyErr_Warn(PyExc_UserWarning,("Warning: getClumpVolume method detected, that clump with id "+lexical_cast<string>(b->getId())+" has non-spherical members.").c_str()); return 0;}
- }
- //generate regular grid points and check if point is inside or outside of clump:
- Vector3r gridPoint;
- int c_in = 0;
- int numPoints1D = ceil(pow(numTries,1./3.));
- for (int ii = 0; ii < numPoints1D; ii++){
- for (int jj = 0; jj < numPoints1D; jj++){
- for (int kk = 0; kk < numPoints1D; kk++){
- bool isInside = false;
- gridPoint = Vector3r(xmin + (xmax-xmin)*(ii+.5)/numPoints1D,ymin + (ymax-ymin)*(jj+.5)/numPoints1D,zmin + (zmax-zmin)*(kk+.5)/numPoints1D);
- FOREACH(MemberMap::value_type& mm, members){
- const Body::id_t& memberId=mm.first;
- const shared_ptr<Body>& member=Body::byId(memberId,scene);
- assert(member->isClumpMember());
- if (member->shape->getClassIndex() == Sph_Index){//clump member should be a sphere
- const Sphere* sphere = YADE_CAST<Sphere*> (member->shape.get());
- if ((gridPoint - member->state->pos).squaredNorm() <= pow(sphere->radius,2.)) { isInside = true; break; }
- }
- }
- if (isInside) c_in += 1;
- }
- }
- }
- //return estimated volume:
- Real vol_estimated = ((1.*c_in)/(1.*pow(numPoints1D,3.)))*abs(xmax-xmin)*abs(ymax-ymin)*abs(zmax-zmin);
- if (vol_estimated <= 0.0) {PyErr_Warn(PyExc_UserWarning,("Something went wrong in getClumpVolume method (clump volume <= 0 detected for clump with id "+lexical_cast<string>(b->getId())+").").c_str()); return 0;}
- else return vol_estimated;
- }
- python::list adaptClumpMasses(python::list excludeList, int num_tries){
- python::list ret;
- Scene* scene(Omega::instance().getScene().get()); // get scene
- shared_ptr<Sphere> sph (new Sphere);
- int Sph_Index = sph->getClassIndexStatic(); // get sphere index for checking if bodies are spheres
- //convert excludeList to a c++ list
- vector<Body::id_t> excludeListC;
- for (int ii = 0; ii < python::len(excludeList); ii++) excludeListC.push_back(python::extract<Body::id_t>(excludeList[ii])());
- Real volClump, estimatedMass, radTmp, volSum, memberMass, inertiaTmp;
- FOREACH(const shared_ptr<Body>& b, *proxee){
- if ( ( !(std::find(excludeListC.begin(), excludeListC.end(), b->getId()) != excludeListC.end()) ) && (b->isClump()) ){
- volClump = getClumpVolume(b,num_tries);
- const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(b->shape);
- std::map<Body::id_t,Se3r>& members = clump->members;
- volSum = 0.0;
- estimatedMass = 0.0;
- for (int ii = 0; ii < 2; ii++){
- FOREACH(MemberMap::value_type& mm, members){
- const Body::id_t& memberId=mm.first;
- const shared_ptr<Body>& member=Body::byId(memberId,scene);
- assert(member->isClumpMember());
- radTmp = 0.0;
- if (member->shape->getClassIndex() == Sph_Index){//clump member should be a sphere
- const Sphere* sphere = YADE_CAST<Sphere*> (member->shape.get());
- radTmp = sphere->radius;
- }
- if (ii == 0) volSum += (4./3.)*Mathr::PI*pow(radTmp,3.);
- else {
- shared_ptr<Material> matTmp = member->material;
- memberMass = (volClump/volSum)*(4./3.)*Mathr::PI*pow(radTmp,3.)*matTmp->density;
- member->state->mass = memberMass;//set vol. corrected mass for clump members
- inertiaTmp = 2.0/5.0*memberMass*radTmp*radTmp;
- member->state->inertia = Vector3r(inertiaTmp,inertiaTmp,inertiaTmp);
- estimatedMass += memberMass;
- }
- }
- }
- Clump::updateProperties(b,/*intersecting*/ false);
- ret.append(python::make_tuple(b->getId(),estimatedMass));
- }
- }
- return ret;
- }
Real getRoundness(python::list excludeList){
Scene* scene(Omega::instance().getScene().get()); // get scene
shared_ptr<Sphere> sph (new Sphere);
@@ -450,13 +351,13 @@
vector<Body::id_t> excludeListC;
for (int ii = 0; ii < python::len(excludeList); ii++) excludeListC.push_back(python::extract<Body::id_t>(excludeList[ii])());
Real RC_sum = 0.0; //sum of local roundnesses
- Real R1, R2, vol_tmp;
+ Real R1, R2, vol, dens;
int c = 0; //counter
FOREACH(const shared_ptr<Body>& b, *proxee){
if ( !(std::find(excludeListC.begin(), excludeListC.end(), b->getId()) != excludeListC.end()) ) {
if ((b->shape->getClassIndex() == Sph_Index) && (b->isStandalone())) { RC_sum += 1.0; c += 1; }
if (b->isClump()){
- R2 = 0.0;
+ R2 = 0.0; dens = 0.0; vol = 0.0;
const shared_ptr<Clump>& clump=YADE_PTR_CAST<Clump>(b->shape);
std::map<Body::id_t,Se3r>& members = clump->members;
FOREACH(MemberMap::value_type& mm, members){
@@ -466,10 +367,11 @@
if (member->shape->getClassIndex() == Sph_Index){//clump member should be a sphere
const Sphere* sphere = YADE_CAST<Sphere*> (member->shape.get());
R2 = max((member->state->pos - b->state->pos).norm() + sphere->radius, R2); //get minimum radius of a sphere, that imbeds clump
+ dens = member->material->density;
}
}
- vol_tmp = getClumpVolume(b,100000);//for RC an accuracy of +/- 1/100000 should be ok
- R1 = pow((3.*vol_tmp)/(4.*Mathr::PI),1./3.); //get theoretical radius of a sphere, with same volume as clump
+ if (dens > 0.) vol = b->state->mass/dens;
+ R1 = pow((3.*vol)/(4.*Mathr::PI),1./3.); //get theoretical radius of a sphere, with same volume as clump
if (R2 < R1) {PyErr_Warn(PyExc_UserWarning,("Something went wrong in getRoundness method (R2 < R1 detected).")); return 0;}
RC_sum += R1/R2; c += 1;
}
@@ -879,9 +781,9 @@
.add_property("dt",&pyOmega::dt_get,&pyOmega::dt_set,"Current timestep (Δt) value.")
.add_property("dynDt",&pyOmega::dynDt_get,&pyOmega::dynDt_set,"Whether a :yref:`TimeStepper` is used for dynamic Δt control. See :yref:`dt<Omega.dt>` on how to enable/disable :yref:`TimeStepper`.")
.add_property("dynDtAvailable",&pyOmega::dynDtAvailable_get,"Whether a :yref:`TimeStepper` is amongst :yref:`O.engines<Omega.engines>`, activated or not.")
- .def("load",&pyOmega::load,(py::arg("file"),py::arg("quiet")=false),"Load simulation from file.")
+ .def("load",&pyOmega::load,(py::arg("file"),py::arg("quiet")=false),"Load simulation from file. The file should be :yref:`saved<Omega.save>` in the same version of Yade, otherwise compatibility is not guaranteed.")
.def("reload",&pyOmega::reload,(py::arg("quiet")=false),"Reload current simulation")
- .def("save",&pyOmega::save,(py::arg("file"),py::arg("quiet")=false),"Save current simulation to file (should be .xml or .xml.bz2)")
+ .def("save",&pyOmega::save,(py::arg("file"),py::arg("quiet")=false),"Save current simulation to file (should be .xml or .xml.bz2). The file should be :yref:`loaded<Omega.load>` in the same version of Yade, otherwise compatibility is not guaranteed.")
.def("loadTmp",&pyOmega::loadTmp,(py::arg("mark")="",py::arg("quiet")=false),"Load simulation previously stored in memory by saveTmp. *mark* optionally distinguishes multiple saved simulations")
.def("saveTmp",&pyOmega::saveTmp,(py::arg("mark")="",py::arg("quiet")=false),"Save simulation to memory (disappears at shutdown), can be loaded later with loadTmp. *mark* optionally distinguishes different memory-saved simulations.")
.def("lsTmp",&pyOmega::lsTmp,"Return list of all memory-saved simulations.")
@@ -938,12 +840,11 @@
.def("__iter__",&pyBodyContainer::pyIter)
.def("append",&pyBodyContainer::append,"Append one Body instance, return its id.")
.def("append",&pyBodyContainer::appendList,"Append list of Body instance, return list of ids")
- .def("appendClumped",&pyBodyContainer::appendClump,"Append given list of bodies as a clump (rigid aggregate); returns a tuple of ``(clumpId,[memberId1,memberId2,...])``.")
- .def("clump",&pyBodyContainer::clump,"Clump given bodies together (creating a rigid aggregate); returns ``clumpId``.")
+ .def("appendClumped",&pyBodyContainer::appendClump,"Append given list of bodies as a clump (rigid aggregate); returns a tuple of ``(clumpId,[memberId1,memberId2,...])``. Clump masses and inertia are adapted automatically. If clump members are overlapping this is done by integration/summation over mass points using a regular grid of cells. For non-overlapping members inertia of the clump is the sum of inertias from members.")
+ .def("clump",&pyBodyContainer::clump,"Clump given bodies together (creating a rigid aggregate); returns ``clumpId``. Clump masses and inertia are adapted automatically (see :yref:`appendClumped()<BodyContainer.appendClumped>`).")
.def("addToClump",&pyBodyContainer::addToClump,"Add body b to an existing clump c. c must be clump and b may not be a clump member of c.\n\nSee **/examples/clumps/addToClump-example.py** for an example script.\n\n.. note:: If b is a clump itself, then all members will be added to c and b will be deleted. If b is a clump member of clump d, then all members from d will be added to c and d will be deleted. If you need to add just clump member b, :yref:`release<BodyContainer.releaseFromClump>` this member from d first.")
.def("releaseFromClump",&pyBodyContainer::releaseFromClump,"Release body b from clump c. b must be a clump member of c.\n\nSee **/examples/clumps/releaseFromClump-example.py** for an example script.\n\n.. note:: If c contains only 2 members b will not be released and a warning will appear. In this case clump c should be :yref:`erased<BodyContainer.erase>`.")
.def("replaceByClumps",&pyBodyContainer::replaceByClumps,"Replace spheres by clumps using a list of clump templates and a list of amounts; returns a list of tuples: ``[(clumpId1,[memberId1,memberId2,...]),(clumpId2,[memberId1,memberId2,...]),...]``. A new clump will have the same volume as the sphere, that was replaced (clump volume/mass/inertia is accounting for overlaps assuming that there are only pair overlaps, to adapt masses of clumps with multiple overlaps use :yref:`adaptClumpMasses()<BodyContainer.adaptClumpMasses>`). \n\n\t *O.bodies.replaceByClumps( [utils.clumpTemplate([1,1],[.5,.5])] , [.9] ) #will replace 90 % of all standalone spheres by 'dyads'*\n\nSee **/examples/clumps/replaceByClumps-example.py** for an example script.")
- .def("adaptClumpMasses",&pyBodyContainer::adaptClumpMasses,"Adapt clump masses and inertia via deterministic Monte-Carlo algorithm, that estimates clump volume using a regular grid of points; returns a list of tuples: ``[(clumpId1,estimatedMassOfClump1),(clumpId2,estimatedMassOfClump2),...]``.\n\nIt is recommended to use this method, when clumps where created via :yref:`clump()<BodyContainer.clump>` or :yref:`appendClumped()<BodyContainer.appendClumped>`. This method could also be used to adapt masses of clumps with multiple overlaps (e.g. triple sections). Bodies can be excluded from the calculation by giving a list of ids: *O.bodies.adaptClumpMasses([ids],1000000) #number of grid points is set to 1000000*\n\nSee **/examples/clumps/adaptClumpMasses-example.py** for an example script.")
.def("getRoundness",&pyBodyContainer::getRoundness,"Returns roundness coefficient RC = R2/R1. R1 is the theoretical radius of a sphere, with same volume as clump. R2 is the minimum radius of a sphere, that imbeds clump. If just spheres are present RC = 1. If clumps are present 0 < RC < 1. Bodies can be excluded from the calculation by giving a list of ids: *O.bodies.getRoundness([ids])*.\n\nSee **/examples/clumps/replaceByClumps-example.py** for an example script.")
.def("clear", &pyBodyContainer::clear,"Remove all bodies (interactions not checked)")
.def("erase", &pyBodyContainer::erase,"Erase body with the given id; all interaction will be deleted by InteractionLoop in the next step.")
=== added file 'scripts/checks-and-tests/checks/checkWirePM.py'
--- scripts/checks-and-tests/checks/checkWirePM.py 1970-01-01 00:00:00 +0000
+++ scripts/checks-and-tests/checks/checkWirePM.py 2013-06-25 03:07:32 +0000
@@ -0,0 +1,145 @@
+# -*- coding: utf-8 -*-
+
+# Check test version for WirePM tensile test
+
+tolerance=0.01
+resultStatus=0
+errors=0
+
+#### define parameters for the net
+# wire diameter
+d = 2.7/1000.
+# particle radius
+radius = d*4.
+# define piecewise lineare stress-strain curve
+strainStressValues=[(0.0019230769,2.5e8),(0.0192,3.2195e8),(0.05,3.8292e8),(0.15,5.1219e8),(0.25,5.5854e8),(0.3,5.6585e8),(0.35,5.6585e8)]
+# elastic material properties
+particleVolume = 4./3.*pow(radius,3)*pi
+particleMass = 3.9/1000.
+density = particleMass/particleVolume
+young = strainStressValues[0][1] / strainStressValues[0][0]
+poisson = 0.3
+
+
+#### material definition
+netMat = O.materials.append( WireMat( young=young,poisson=poisson,frictionAngle=radians(30),density=density,isDoubleTwist=True,diameter=d,strainStressValues=strainStressValues,lambdaEps=0.4,lambdak=0.66) )
+
+wireMat = O.materials.append( WireMat( young=young,poisson=poisson,frictionAngle=radians(30),density=density,isDoubleTwist=False,diameter=3.4/1000,strainStressValues=strainStressValues ) )
+
+
+#### get net packing
+kw = {'color':[1,1,0],'wire':True,'highlight':False,'fixed':False,'material':netMat}
+[netpack,lx,ly] = hexaNet( radius=radius, cornerCoord=[0,0,0], xLength=1.0, yLength=0.55, mos=0.08, a=0.04, b=0.04, startAtCorner=False, isSymmetric=True, **kw )
+O.bodies.append(netpack)
+
+if abs((lx-0.96)/0.96)>tolerance:
+ print "WirePM checkTest: difference on net dimension lx"
+ errors+=1
+if abs((ly-0.52)/0.52)>tolerance:
+ print "WirePM checkTest: difference on net dimension ly"
+ errors+=1
+
+
+#### get bodies for single wire at the boundary in y-direction and change properties
+bb = uniaxialTestFeatures(axis=0)
+negIds,posIds=bb['negIds'],bb['posIds']
+
+for id in negIds:
+ O.bodies[id].material = O.materials[wireMat]
+ O.bodies[id].shape.color = [0,0,1]
+for id in posIds:
+ O.bodies[id].material = O.materials[wireMat]
+ O.bodies[id].shape.color = [0,0,1]
+
+
+#### define engines to create link
+interactionRadius=2.8
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='aabb')]),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='Ig2ssGeom')],
+ [Ip2_WireMat_WireMat_WirePhys(linkThresholdIteration=1,label='interactionPhys')],
+ [Law2_ScGeom_WirePhys_WirePM(linkThresholdIteration=1,label='interactionLaw')]
+ ),
+ NewtonIntegrator(damping=0.),
+]
+
+#### define additional vertical interactions at the boundary
+createInteraction(negIds[0],negIds[2])
+createInteraction(negIds[3],negIds[4])
+createInteraction(negIds[5],negIds[6])
+createInteraction(negIds[7],negIds[1])
+createInteraction(posIds[0],posIds[2])
+createInteraction(posIds[3],posIds[4])
+createInteraction(posIds[5],posIds[6])
+createInteraction(posIds[7],posIds[1])
+
+#### time step definition for first time step to create links
+O.step()
+
+
+#### initialize values for UniaxialStrainer
+bb = uniaxialTestFeatures(axis=1)
+negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
+strainRateTension = 0.1
+setSpeeds = True
+
+
+##### delete horizontal interactions for corner particles
+bb = uniaxialTestFeatures(axis=1)
+negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
+
+
+##### delete some interactions
+O.interactions.erase(0,4)
+O.interactions.erase(0,5)
+O.interactions.erase(1,154)
+O.interactions.erase(1,155)
+O.interactions.erase(2,26)
+O.interactions.erase(2,27)
+O.interactions.erase(3,176)
+O.interactions.erase(3,177)
+
+#### time step definition for deleting some links which have been created by the Ig2 functor
+O.step()
+
+
+#### initializes now the interaction detection factor
+aabb.aabbEnlargeFactor=-1.
+Ig2ssGeom.interactionDetectionFactor=-1.
+
+
+#### define engines for simulation with UniaxialStrainer
+O.engines = O.engines[:3] + [ UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=1,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=True,blockRotations=False,setSpeeds=setSpeeds,label='strainer'),
+NewtonIntegrator(damping=0.5),
+]
+
+
+#### time step definition for simulation
+## critical time step proposed by Bertrand
+kn = 16115042 # stiffness of single wire from code
+O.dt = 0.2*sqrt(particleMass/(2.*kn))
+
+O.run(350000,True)
+
+Fn = 0.
+for i in posIds:
+ try:
+ inter=O.interactions.withBody(i)[0]
+ F = abs(inter.phys.normalForce[1])
+ except:
+ F = 0
+ Fn += F
+un = O.bodies[O.bodies[posIds[0]].id].state.pos[1] - O.bodies[O.bodies[posIds[0]].id].state.refPos[1]
+
+
+if abs((un-0.04)/0.04)>tolerance :
+ print "WirePM checkTest: difference on peak displacement"
+ errors+=1
+if abs((Fn-33626.44)/33626.44)>tolerance :
+ print "WirePM checkTest: difference on peak Force"
+ errors+=1
+
+if (errors):
+ resultStatus +=1 #Test is failed
