TetGen and third party libraries

Hi Hicham,

there is a patch, which is linking TETGEN:

Subject: [PATCH] Add TetGen coupling.

---
 CMakeLists.txt | 17 +++++++++
 cMake/FindTetGen.cmake | 17 +++++++++
 pkg/dem/Shop.hpp | 8 ++++
 pkg/dem/Shop_02.cpp | 102 +++++++++++++++++++++++++++++++++++++++++++++++++
 py/_utils.cpp | 2 +
 5 files changed, 146 insertions(+)
 create mode 100644 cMake/FindTetGen.cmake

diff --git a/CMakeLists.txt b/CMakeLists.txt
index a299e46..70b4a92 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,6 +18,7 @@
 # ENABLE_LBMFLOW: enable LBMFLOW-option, LBM_ENGINE (OFF by default)
 # ENABLE_LIQMIGRATION: enable LIQMIGRATION-option, see [Mani2013]
for details (OFF by default)
 # ENABLE_MASK_ARBITRARY: enable MASK_ARBITRARY option (OFF by default)
+# ENABLE_TETGEN: enable TETGEN option (OFF by default)
 # runtimePREFIX: used for packaging, when install directory is not
the same is runtime directory (/usr/local by default)
 # CHUNKSIZE: set >1, if you want several sources to be compiled at
once. Increases compilation speed and RAM-consumption during it (1 by
default).
 # VECTORIZE: enables vectorization and alignment in Eigen3 library,
experimental (OFF by default)
@@ -126,6 +127,7 @@ OPTION(ENABLE_PFVFLOW "Enable flow engine
(experimental)" ${DEFAULT})
 OPTION(ENABLE_SPH "Enable SPH" OFF)
 OPTION(ENABLE_LBMFLOW "Enable LBM engine (very experimental)" OFF)
 OPTION(ENABLE_LIQMIGRATION "Enable liquid control (very
experimental), see [Mani2013] for details" OFF)
+OPTION(ENABLE_TETGEN "Enable TetGen" OFF)
 OPTION(ENABLE_MASK_ARBITRARY "Enable arbitrary precision of bitmask
variables (only Body::groupMask yet implemented) (experimental). Use
-DMASK_ARBITRARY_SIZE=int to set number of used bits (256 by default)"
OFF)

 #===========================================================
@@ -320,6 +322,16 @@ ELSE(ENABLE_LIQMIGRATION)
   SET(DISABLED_FEATS "${DISABLED_FEATS} LIQMIGRATION")
 ENDIF(ENABLE_LIQMIGRATION)
 #===============================================
+IF(ENABLE_TETGEN)
+ FIND_PACKAGE(TetGen REQUIRED)
+ ADD_DEFINITIONS("-DTETGEN")
+ ADD_DEFINITIONS(${TETGEN_DEFINITIONS})
+ INCLUDE_DIRECTORIES(${TETGEN_INCLUDE_DIR})
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} TETGEN")
+ELSE(ENABLE_TETGEN)
+ SET(DISABLED_FEATS "${DISABLED_FEATS} TETGEN")
+ENDIF(ENABLE_TETGEN)
+#===============================================
 IF(ENABLE_GL2PS)
   FIND_PACKAGE(GL2PS)
   IF(GL2PS_FOUND)
@@ -478,6 +490,11 @@ IF(ENABLE_GUI)
 ENDIF(ENABLE_GUI)

 #====================================
+IF(ENABLE_TETGEN)
+ TARGET_LINK_LIBRARIES(yade ${TETGEN_LIBRARIES})
+ENDIF(ENABLE_TETGEN)
+#====================================
+
 IF (NOT (PY_minieigen))
   ADD_LIBRARY(miniEigen SHARED py/mathWrap/miniEigen.cpp)
   SET_TARGET_PROPERTIES(miniEigen PROPERTIES PREFIX "")
diff --git a/cMake/FindTetGen.cmake b/cMake/FindTetGen.cmake
new file mode 100644
index 0000000..ab8c0a6
--- /dev/null
+++ b/cMake/FindTetGen.cmake
@@ -0,0 +1,17 @@
+# Try to find the TetGen librairies
+# TETGEN_FOUND - system has TetGen lib
+# TETGEN_INCLUDE_DIR - the TetGen include directory
+# TETGEN_LIBRARIES - Libraries needed to use TetGen
+# TETGEN_DEFINITIONS - Compiler switches required for using gnutls
+
+IF (TETGEN_INCLUDE_DIR AND TETGEN_LIBRARIES)
+ SET(TETGEN_FIND_QUIETLY TRUE)
+ENDIF (TETGEN_INCLUDE_DIR AND TETGEN_LIBRARIES)
+
+SET(TETGEN_DEFINITIONS -DTETLIBRARY)
+FIND_PATH(TETGEN_INCLUDE_DIR NAMES tetgen.h )
+FIND_LIBRARY(TETGEN_LIBRARIES NAMES tet )
+MESSAGE(STATUS "TetGen lib: " ${TETGEN_LIBRARIES})
+
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(TETGEN DEFAULT_MSG
TETGEN_INCLUDE_DIR TETGEN_LIBRARIES)
diff --git a/pkg/dem/Shop.hpp b/pkg/dem/Shop.hpp
index ef45538..f766971 100644
--- a/pkg/dem/Shop.hpp
+++ b/pkg/dem/Shop.hpp
@@ -15,6 +15,9 @@

 #include<boost/function.hpp>

+#ifdef TETGEN
+#include "tetgen.h"
+#endif

 class Scene;
 class Body;
@@ -148,4 +151,9 @@ class Shop{
         static void growParticles(Real multiplier, bool updateMass,
bool dynamicOnly, unsigned int discretization);
         //! Change of size of a single sphere or a clump
         static void growParticle(Body::id_t bodyID, Real multiplier,
bool updateMass);
+
+#ifdef TETGEN
+ //! Test function to show how TetGen works
+ static py::list addTetgenBody();
+#endif
 };
diff --git a/pkg/dem/Shop_02.cpp b/pkg/dem/Shop_02.cpp
index 3182128..02dc616 100644
--- a/pkg/dem/Shop_02.cpp
+++ b/pkg/dem/Shop_02.cpp
@@ -53,6 +53,11 @@
     #define FOREACH BOOST_FOREACH
 #endif

+#ifdef TETGEN
+#include <yade/pkg/common/Facet.hpp>
+#include <yade/core/Scene.hpp>
+#endif
+
 CREATE_LOGGER(Shop);

 /*! Flip periodic cell by given number of cells.
@@ -515,3 +520,100 @@ void Shop::growParticle(Body::id_t bodyID, Real
multiplier, bool updateMass)
         else contact->refR2 = rad;
     }
 }
+
+py::list Shop::addTetgenBody()
+{
+ const Scene* scene=Omega::instance().getScene().get();
+ py::list ret;
+ tetgenio in, out;
+ tetgenio::facet *f;
+ tetgenio::polygon *p;
+
+ // All indices start from 1.
+ in.firstnumber = 1;
+
+ // Points
+ std::vector<Eigen::Vector3d> pointsV;
+
+ pointsV.push_back(Eigen::Vector3d(0,0,0)); // node 1.
+ pointsV.push_back(Eigen::Vector3d(2,0,0)); // node 2.
+ pointsV.push_back(Eigen::Vector3d(2,2,0)); // node 3.
+ pointsV.push_back(Eigen::Vector3d(0,2,0)); // node 4.
+ pointsV.push_back(Eigen::Vector3d(1,1,5)); // node 5.
+
+ // Facets
+ std::vector<Eigen::Vector3i> facetsV;
+ facetsV.push_back(Eigen::Vector3i(1,2,3)); // facet 1.
+ facetsV.push_back(Eigen::Vector3i(1,3,4)); // facet 2.
+ facetsV.push_back(Eigen::Vector3i(1,2,5)); // facet 3.
+ facetsV.push_back(Eigen::Vector3i(2,3,5)); // facet 4.
+ facetsV.push_back(Eigen::Vector3i(3,4,5)); // facet 5.
+ facetsV.push_back(Eigen::Vector3i(4,5,1)); // facet 6.
+
+ // Add points
+ in.numberofpoints = pointsV.size();
+ in.pointlist = new REAL[in.numberofpoints * 3];
+ for (int d = 0; d < in.numberofpoints; d++) {
+ in.pointlist[3*d+0] = pointsV[d][0];
+ in.pointlist[3*d+1] = pointsV[d][1];
+ in.pointlist[3*d+2] = pointsV[d][2];
+ }
+
+ // Add facets
+ in.numberoffacets = facetsV.size();
+ in.facetlist = new tetgenio::facet[in.numberoffacets];
+ in.facetmarkerlist = new int[in.numberoffacets];
+
+ for (int d = 0; d < in.numberoffacets; d++) {
+ f = &in.facetlist[d];
+ f->numberofpolygons = 1;
+ f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
+ f->numberofholes = 0;
+ f->holelist = NULL;
+ p = &f->polygonlist[0];
+ p->numberofvertices = 3;
+ p->vertexlist = new int[p->numberofvertices];
+ p->vertexlist[0] = facetsV[d][0];
+ p->vertexlist[1] = facetsV[d][1];
+ p->vertexlist[2] = facetsV[d][2];
+ }
+
+ tetrahedralize("pq1.4a0.1", &in, &out);
+
+ std::list<Vector3i> allTriag;
+ for (int d = 0; d < out.numberoftetrahedra; d++) {
+ std::vector<int> v;
+ v.push_back(out.tetrahedronlist[d*4 + 0]);
+ v.push_back(out.tetrahedronlist[d*4 + 1]);
+ v.push_back(out.tetrahedronlist[d*4 + 2]);
+ v.push_back(out.tetrahedronlist[d*4 + 3]);
+ std::sort (v.begin(), v.end());
+ allTriag.push_back(Vector3i(v[0]-1,v[1]-1,v[2]-1));
+ allTriag.push_back(Vector3i(v[0]-1,v[1]-1,v[3]-1));
+ allTriag.push_back(Vector3i(v[0]-1,v[2]-1,v[3]-1));
+ allTriag.push_back(Vector3i(v[1]-1,v[2]-1,v[3]-1));
+ }
+ allTriag.unique();
+ std::cout<<"Tiangles :"<<allTriag.size()<<std::endl;
+ for (std::list<Vector3i>::iterator it=allTriag.begin(); it !=
allTriag.end(); ++it) {
+ Vector3r v[3];
+ v[0]=Vector3r(out.pointlist[(*it)[0]*3 +
0],out.pointlist[(*it)[0]*3 + 1],out.pointlist[(*it)[0]*3 + 2]);
+ v[1]=Vector3r(out.pointlist[(*it)[1]*3 +
0],out.pointlist[(*it)[1]*3 + 1],out.pointlist[(*it)[1]*3 + 2]);
+ v[2]=Vector3r(out.pointlist[(*it)[2]*3 +
0],out.pointlist[(*it)[2]*3 + 1],out.pointlist[(*it)[2]*3 + 2]);
+
+ Vector3r icc = Shop::inscribedCircleCenter(v[0],v[1],v[2]);
+ shared_ptr<Facet> iFacet(new Facet);
+ iFacet->color = Vector3r(0.8,0.3,0.3);
+ for (int j=0; j<3; ++j) {
+ iFacet->vertices[j]=v[j]-icc;
+ }
+ iFacet->postLoad(*iFacet);
+ shared_ptr<Body> b(new Body);
+ b->state->pos=b->state->refPos=icc;
+ b->state->ori=b->state->refOri=Quaternionr::Identity();
+ b->shape = iFacet;
+ b->shape->wire = true;
+ ret.append(scene->bodies->insert(b));
+ }
+ return ret;
+}
diff --git a/py/_utils.cpp b/py/_utils.cpp
index 2de0a43..019f815 100644
--- a/py/_utils.cpp
+++ b/py/_utils.cpp
@@ -389,6 +389,7 @@ py::tuple Shop__normalShearStressTensors(bool
compressionPositive=false, bool sp
 py::list Shop__getStressLWForEachBody(){return Shop::getStressLWForEachBody();}

 py::list Shop__getBodyIdsContacts(Body::id_t bodyID=-1){return
Shop::getBodyIdsContacts(bodyID);}
+py::list Shop__addTetgenBody(){return Shop::addTetgenBody();}

 Real shiftBodies(py::list ids, const Vector3r& shift){
     shared_ptr<Scene> rb=Omega::instance().getScene();
@@ -495,6 +496,7 @@ BOOST_PYTHON_MODULE(_utils){
     py::def("getStress",Shop::getStress,(py::args("volume")=0),"Compute
and return Love-Weber stress tensor:\n\n
$\\sigma_{ij}=\\frac{1}{V}\\sum_b f_i^b l_j^b$, where the sum is over
all interactions, with $f$ the contact force and $l$ the branch vector
(joining centers of the bodies). Stress is negativ for repulsive
contact forces, i.e. compression. $V$ can be passed to the function.
If it is not, it will be equal to one in non-periodic cases, or equal
to the volume of the cell in periodic cases.");
     py::def("getCapillaryStress",Shop::getCapillaryStress,(py::args("volume")=0),"Compute
and return Love-Weber capillary stress tensor:\n\n
$\\sigma^{cap}_{ij}=\\frac{1}{V}\\sum_b l_i^b f^{cap,b}_j$, where the
sum is over all interactions, with $l$ the branch vector (joining
centers of the bodies) and $f^{cap}$ is the capillary force. $V$ can
be passed to the function. If it is not, it will be equal to one in
non-periodic cases, or equal to the volume of the cell in periodic
cases. Only the CapillaryPhys interaction type is supported
presently.");
     py::def("getBodyIdsContacts",Shop__getBodyIdsContacts,(py::args("bodyID")=0),"Get
a list of body-ids, which contacts the given body.");
+ py::def("addTetgenBody",Shop__addTetgenBody,"Add TetGen body.");
     py::def("maxOverlapRatio",maxOverlapRatio,"Return maximum overlap
ration in interactions (with :yref:`ScGeom`) of two
:yref:`spheres<Sphere>`. The ratio is computed as $\\frac{u_N}{2(r_1
r_2)/r_1+r_2}$, where $u_N$ is the current overlap distance and $r_1$,
$r_2$ are radii of the two spheres in contact.");
     py::def("shiftBodies",shiftBodies,(py::arg("ids"),py::arg("shift")),"Shifts
bodies listed in ids without updating their velocities.");
     py::def("calm",Shop__calm,(py::arg("mask")=-1),"Set translational
and rotational velocities of bodies to zero. Applied to all bodies by
default. To calm only some bodies, use mask parameter, it will calm
only bodies with groupMask compatible to given value");
--
2.0.0

Anton

2014-07-10 13:36 GMT+02:00 Jan Stránský <email address hidden>:
> Question #251458 on Yade changed:
> https://answers.launchpad.net/yade/+question/251458
>
> Status: Open => Answered
>
> Jan Stránský proposed the following answer:
> Hi Hicham,
> thanks for the info. Did you define the tetrahedralize function in some
> .cpp file inside BOOST_PYTHON_MODULE, and if so, where and how?
> thanks
> Jan
>
>
>
> 2014-07-10 13:17 GMT+02:00 Hicham BENNIOU <
> <email address hidden>>:
>
>> Question #251458 on Yade changed:
>> https://answers.launchpad.net/yade/+question/251458
>>
>> Status: Answered => Open
>>
>> Hicham BENNIOU is still having a problem:
>> Hi Jan,
>>
>> I get no compilation error, it goes smoothly and YADE runs correctly,
>> but when I try to call "tetrahedralize()" function for example I get an
>> error :
>>
>> NameError: name 'tetrahedralize' is not defined
>>
>> --
>> You received this question notification because you are a member of
>> yade-users, which is an answer contact for Yade.
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~yade-users
>> Post to : <email address hidden>
>> Unsubscribe : https://launchpad.net/~yade-users
>> More help : https://help.launchpad.net/ListHelp
>>
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp