yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12483
[Branch ~yade-pkg/yade/git-trunk] Rev 3760: Initial version of the PFacet implementation (contributed by Anna Effeindzourou)
------------------------------------------------------------
revno: 3760
committer: Klaus Thoeni <klaus.thoeni@xxxxxxxxx>
timestamp: Tue 2015-12-15 23:57:49 +1100
message:
Initial version of the PFacet implementation (contributed by Anna Effeindzourou)
added:
pkg/common/Gl1_PFacet.cpp
pkg/common/Gl1_PFacet.hpp
pkg/common/PFacet.cpp
pkg/common/PFacet.hpp
modified:
pkg/common/Grid.cpp
pkg/common/Grid.hpp
py/gridpfacet.py
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== added file 'pkg/common/Gl1_PFacet.cpp'
--- pkg/common/Gl1_PFacet.cpp 1970-01-01 00:00:00 +0000
+++ pkg/common/Gl1_PFacet.cpp 2015-12-15 12:57:49 +0000
@@ -0,0 +1,57 @@
+#ifdef YADE_OPENGL
+#include<pkg/common/Gl1_PFacet.hpp>
+
+#include<lib/opengl/OpenGLWrapper.hpp>
+
+// bool Gl1_PFacet::normals=false;
+bool Gl1_PFacet::wire=true;
+void Gl1_PFacet::go(const shared_ptr<Shape>& cm, const shared_ptr<State>& ,bool wire2,const GLViewInfo&)
+{
+ PFacet* Pfacet = static_cast<PFacet*>(cm.get());
+ vector<Vector3r> vertices;
+ vertices.push_back(Pfacet->node1->state->pos);
+ vertices.push_back(Pfacet->node2->state->pos);
+ vertices.push_back(Pfacet->node3->state->pos);
+
+ Vector3r pos=Pfacet->node1->state->pos;
+
+ vertices[0]=vertices[0]-pos;
+ vertices[1]=vertices[1]-pos;
+ vertices[2]=vertices[2]-pos;
+
+ vector<Vector3r> verticesF1 = vertices;
+ Vector3r normal=(vertices[1]-vertices[0]).cross(vertices[2]-vertices[1]); normal.normalize();
+ verticesF1[0]=vertices[0] + normal*Pfacet->radius;
+ verticesF1[1]=vertices[1] + normal*Pfacet->radius;
+ verticesF1[2]=vertices[2] + normal*Pfacet->radius;
+
+ vector<Vector3r> verticesF2 = vertices;
+
+ verticesF2[0] = vertices[0] - normal*Pfacet->radius;
+ verticesF2[1] = vertices[1] - normal*Pfacet->radius;
+ verticesF2[2] = vertices[2] - normal*Pfacet->radius;
+
+ if(!wire2||!wire){
+
+ glDisable(GL_CULL_FACE);
+
+ glColor3v(cm->color);
+ glBegin(GL_TRIANGLES);
+ glNormal3v(normal); // this makes every triangle different WRT the light direction; important!
+ glVertex3v(verticesF1[0]);
+ glVertex3v(verticesF1[1]);
+ glVertex3v(verticesF1[2]);
+ glEnd();
+ glBegin(GL_TRIANGLES);
+ glNormal3v(Pfacet->normal); // this makes every triangle different WRT the light direction; important!
+ glVertex3v(verticesF2[2]);
+ glVertex3v(verticesF2[1]);
+ glVertex3v(verticesF2[0]);
+ glEnd();
+
+ }
+}
+
+YADE_PLUGIN((Gl1_PFacet));
+
+#endif
\ No newline at end of file
=== added file 'pkg/common/Gl1_PFacet.hpp'
--- pkg/common/Gl1_PFacet.hpp 1970-01-01 00:00:00 +0000
+++ pkg/common/Gl1_PFacet.hpp 2015-12-15 12:57:49 +0000
@@ -0,0 +1,21 @@
+#pragma once
+#include<pkg/common/PFacet.hpp>
+#include<pkg/common/GLDrawFunctors.hpp>
+#include<pkg/common/Facet.hpp>
+#include<core/Shape.hpp>
+
+#ifdef YADE_OPENGL
+class Gl1_PFacet : public GlShapeFunctor
+{
+ public:
+ virtual void go(const shared_ptr<Shape>&, const shared_ptr<State>&,bool,const GLViewInfo&);
+ RENDERS(PFacet);
+ YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_PFacet,GlShapeFunctor,"Renders :yref:`Facet` object",
+ ((bool,wire,false,,"Only show wireframe (controlled by ``glutSlices`` and ``glutStacks``."))
+ );
+};
+
+REGISTER_SERIALIZABLE(Gl1_PFacet);
+#endif
+
+
=== modified file 'pkg/common/Grid.cpp'
--- pkg/common/Grid.cpp 2014-09-23 13:00:45 +0000
+++ pkg/common/Grid.cpp 2015-12-15 12:57:49 +0000
@@ -39,7 +39,46 @@
return getSegment().norm();
}
-
+void GridNode::addPFacet(shared_ptr<Body> PF){
+ pfacetList.push_back(PF);
+}
+
+void GridConnection::addPFacet(shared_ptr<Body> PF){
+ pfacetList.push_back(PF);
+}
+
+PFacet::~PFacet(){}
+YADE_PLUGIN((PFacet));
+void PFacet::postLoad(PFacet&)
+{
+ vector<Vector3r> vertices;
+ vertices.push_back(node1->state->pos);
+ vertices.push_back(node2->state->pos);
+ vertices.push_back(node3->state->pos);
+
+ Vector3r center;
+ center=vertices[0]+((vertices[2]-vertices[0])*(vertices[1]-vertices[0]).norm()+(vertices[1]-vertices[0])*(vertices[2]-vertices[0]).norm())/((vertices[1]-vertices[0]).norm()+(vertices[2]-vertices[1]).norm()+(vertices[0]-vertices[2]).norm());
+
+ vertices[0]=vertices[0]-center;
+ vertices[1]=vertices[1]-center;
+ vertices[2]=vertices[2]-center;
+ if(vertices.size()!=3){ throw runtime_error(("Facet must have exactly 3 vertices (not "+boost::lexical_cast<string>(vertices.size())+")").c_str()); }
+ if(isnan(vertices[0][0])) return; // not initialized, nothing to do
+ Vector3r e[3] = {vertices[1]-vertices[0] ,vertices[2]-vertices[1] ,vertices[0]-vertices[2]};
+ #define CHECK_EDGE(i) if(e[i].squaredNorm()==0){LOG_FATAL("Facet has coincident vertices "<<i<<" ("<<vertices[i]<<") and "<<(i+1)%3<<" ("<<vertices[(i+1)%3]<<")!");}
+ CHECK_EDGE(0); CHECK_EDGE(1);CHECK_EDGE(2);
+ #undef CHECK_EDGE
+ normal = e[0].cross(e[1]);
+ area = .5*normal.norm();
+ normal /= 2*area;
+ for(int i=0; i<3; ++i){
+ ne[i]=e[i].cross(normal); ne[i].normalize();
+ vl[i]=vertices[i].norm();
+ vu[i]=vertices[i]/vl[i];
+ }
+ Real p = e[0].norm()+e[1].norm()+e[2].norm();
+ icr = e[0].norm()*ne[0].dot(e[2])/p;
+}
//!################## IGeom Functors #####################
@@ -48,6 +87,7 @@
{
//GridConnection* GC = static_cast<GridConnection*>(cm.get());
bool isNew = !c->geom;
+ GridNode* GN[2]={static_cast<GridNode*>(cm1.get()),static_cast<GridNode*>(cm2.get())};
if (Ig2_Sphere_Sphere_ScGeom::go(cm1,cm2,state1,state2,shift2,force,c)){//the 3 DOFS from ScGeom are updated here
if (isNew) {//generate a 6DOF interaction from the 3DOF one generated by Ig2_Sphere_Sphere_ScGeom
shared_ptr<GridNodeGeom6D> sc (new GridNodeGeom6D());
@@ -57,6 +97,14 @@
if (updateRotations) YADE_PTR_CAST<GridNodeGeom6D>(c->geom)->precomputeRotations(state1,state2,isNew,creep);
if(YADE_PTR_CAST<GridNodeGeom6D>(c->geom)->connectionBody){ //test this because the connectionBody may not have been yet initialized.
YADE_PTR_CAST<GridNodeGeom6D>(c->geom)->connectionBody->state->pos=state1.pos;
+
+ for (unsigned int j=0; j<2; j++){
+ for (unsigned int i=0; i< GN[j]->pfacetList.size() ;i++){
+ PFacet* Pfacet = YADE_CAST<PFacet*>(GN[j]->pfacetList[i]->shape.get());
+ if(c->id1==Pfacet->node1->getId())
+ GN[j]->pfacetList[i]->state->pos=state1.pos;
+ }
+ }
}
return true;
}
@@ -421,12 +469,15 @@
}
Real& un=geom->penetrationDepth;
phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
-
+// std::cout<< "phys->shearForce1= "<<phys->shearForce<<std::endl;
Vector3r& shearForce = geom->rotate(phys->shearForce);
+// std::cout<< "phys->shearForce2= "<<phys->shearForce<<", geom->shearIncrement()= "<<geom->shearIncrement()<<", phys->ks= "<<phys->ks<<std::endl;
const Vector3r& shearDisp = geom->shearIncrement();
shearForce -= phys->ks*shearDisp;
Real maxFs = phys->normalForce.squaredNorm()*std::pow(phys->tangensOfFrictionAngle,2);
-
+
+// std::cout<< "shearForce3= "<<shearForce<<" ,phys->normalForce="<<phys->normalForce<<std::endl;
+
if (!scene->trackEnergy){//Update force but don't compute energy terms (see below))
// PFC3d SlipModel, is using friction angle. CoulombCriterion
if( shearForce.squaredNorm() > maxFs ){
@@ -446,13 +497,24 @@
scene->energy->add(0.5*(phys->normalForce.squaredNorm()/phys->kn+phys->shearForce.squaredNorm()/phys->ks),"elastPotential",elastPotentialIx,/*reset at every timestep*/true);
}
Vector3r force = -phys->normalForce-shearForce;
+// std::cout<< "id1= "<<id1<<" ,force="<<force<<std::endl;
scene->forces.addForce(id1,force);
scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
Vector3r twist = (geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force);
- scene->forces.addForce(geom->id3,(geom->relPos-1)*force);
- scene->forces.addTorque(geom->id3,(1-geom->relPos)*twist);
- scene->forces.addForce(geom->id4,(-geom->relPos)*force);
- scene->forces.addTorque(geom->id4,geom->relPos*twist);
+
+
+ if (geom->id5==-1){
+ scene->forces.addForce(geom->id3,(geom->relPos-1)*force);
+ scene->forces.addTorque(geom->id3,(1-geom->relPos)*twist);
+ scene->forces.addForce(geom->id4,(-geom->relPos)*force);
+ scene->forces.addTorque(geom->id4,geom->relPos*twist);
+ }
+ else{
+ scene->forces.addForce(geom->id3,geom->weight[0]*-force);
+ scene->forces.addForce(geom->id4,geom->weight[1]*-force);
+ scene->forces.addForce(geom->id5,geom->weight[2]*-force);
+ }
+
return true;
}
YADE_PLUGIN((Law2_ScGridCoGeom_FrictPhys_CundallStrack));
=== modified file 'pkg/common/Grid.hpp'
--- pkg/common/Grid.hpp 2014-10-15 06:44:01 +0000
+++ pkg/common/Grid.hpp 2015-12-15 12:57:49 +0000
@@ -18,8 +18,6 @@
- The Law2_ScGridCoGeom_FrictPhys_CundallStrack who handles the elastic frictional Sphere-GridConnection contact. The GridNode-GridNode law is Law2_ScGeom6D_CohFrictPhys_CohesionMoment by inheritance.
*/
-
-
#pragma once
#include "Sphere.hpp"
#include <pkg/dem/FrictPhys.hpp>
@@ -42,12 +40,16 @@
virtual ~GridConnection();
Real getLength();
Vector3r getSegment();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(GridConnection,Sphere,"GridConnection shape. Component of a grid designed to link two :yref:`GridNodes<GridNode>`. It's highly recommended to use utils.gridConnection(...) to generate correct :yref:`GridConnections<GridConnection>`.",
+ void addPFacet(shared_ptr<Body> PF);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(GridConnection,Sphere,"GridConnection shape. Component of a grid designed to link two :yref:`GridNodes<GridNode>`. It's highly recommended to use utils.gridConnection(...) to generate correct :yref:`GridConnections<GridConnection>`.",
((shared_ptr<Body> , node1 , ,,"First :yref:`Body` the GridConnection is connected to."))
((shared_ptr<Body> , node2 , ,,"Second :yref:`Body` the GridConnection is connected to."))
((bool, periodic, false,,"true if two nodes from different periods are connected."))
+ ((vector<shared_ptr<Body> >,pfacetList,,,"List of :yref:`PFacet<PFacet>` the GridConnection is connected to."))
((Vector3i , cellDist , Vector3i(0,0,0),,"missing doc :(")),
- createIndex(); /*ctor*/
+ createIndex();, /*ctor*/
+ /*py*/
+ .def("addPFacet",&GridConnection::addPFacet,(boost::python::arg("Body")),"Add a PFacet to the GridConnection.")
);
REGISTER_CLASS_INDEX(GridConnection,Sphere);
};
@@ -58,18 +60,58 @@
public:
virtual ~GridNode();
void addConnection(shared_ptr<Body> GC);
+ void addPFacet(shared_ptr<Body> PF);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(GridNode,Sphere,"GridNode shape, component of a grid.\nTo create a Grid, place the nodes first, they will define the spacial discretisation of it. It's highly recommended to use utils.gridNode(...) to generate correct :yref:`GridNodes<GridNode>`. Note that the GridNodes should only be in an Interaction with other GridNodes. The Sphere-Grid contact is only handled by the :yref:`GridConnections<GridConnection>`.",
+ ((vector<shared_ptr<Body> >,pfacetList,,,"List of :yref:`PFacet<PFacet>` the GridConnection is connected to."))
((vector<shared_ptr<Body> >,ConnList,,,"List of :yref:`GridConnections<GridConnection>` the GridNode is connected to.")),
/*ctor*/
createIndex();,
/*py*/
.def("addConnection",&GridNode::addConnection,(boost::python::arg("Body")),"Add a GridConnection to the GridNode.")
+ .def("addPFacet",&GridNode::addPFacet,(boost::python::arg("Body")),"Add a PFacet to the GridNode.")
);
REGISTER_CLASS_INDEX(GridNode,Sphere);
};
REGISTER_SERIALIZABLE(GridNode);
+//!################## PFacet SHAPES #####################
+class PFacet : public Shape {
+ public:
+
+ virtual ~PFacet();
+ /// Normals of edges
+ Vector3r ne[3];
+ /// Inscribing cirle radius
+ Real icr;
+ /// Length of the vertice vectors
+ Real vl[3];
+ /// Unit vertice vectors
+ Vector3r vu[3];
+ void postLoad(PFacet&);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(PFacet,Shape,"PFacet (particle facet) geometry.",
+ ((shared_ptr<Body> , node1 , ,,"First :yref:`Body` the Pfacet is connected to."))
+ ((shared_ptr<Body> , node2 , ,,"Second :yref:`Body` the Pfacet is connected to."))
+ ((shared_ptr<Body> , node3 , ,,"third :yref:`Body` the Pfacet is connected to."))
+ ((shared_ptr<Body> , conn1 , ,,"First :yref:`Body` the Pfacet is connected to."))
+ ((shared_ptr<Body> , conn2 , ,,"Second :yref:`Body` the Pfacet is connected to."))
+ ((shared_ptr<Body> , conn3 , ,,"third :yref:`Body` the Pfacet is connected to."))
+ ((Vector3r,normal,Vector3r(NaN,NaN,NaN),(Attr::readonly | Attr::noSave),"PFacet's normal (in local coordinate system)"))
+ ((Real,radius,-1,,"PFacet's radius"))
+ ((Real,area,NaN,(Attr::readonly | Attr::noSave),"PFacet's area"))
+ ((Vector3i , cellDist , Vector3i(0,0,0),,"missing doc :("))
+ #ifdef FACET_TOPO
+ ((vector<Body::id_t>,edgeAdjIds,vector<Body::id_t>(3,Body::ID_NONE),,"PFacet id's that are adjacent to respective edges [experimental]"))
+ ((vector<Real>,edgeAdjHalfAngle,vector<Real>(3,0),,"half angle between normals of this facet and the adjacent facet [experimental]"))
+ #endif
+ ,
+ /* ctor */ createIndex();
+ );
+ DECLARE_LOGGER;
+
+ REGISTER_CLASS_INDEX(PFacet,Shape);
+};
+REGISTER_SERIALIZABLE(PFacet);
//!################## Contact Geometry #####################
//! O-O
@@ -97,6 +139,8 @@
((int,trueInt,-1,,"Defines the body id of the :yref:`GridConnection` where the contact is real, when :yref:`ScGridCoGeom::isDuplicate`>0."))
((int,id3,0,,"id of the first :yref:`GridNode`. |yupdate|"))
((int,id4,0,,"id of the second :yref:`GridNode`. |yupdate|"))
+ ((int,id5,-1,,"id of the third :yref:`GridNode`. |yupdate|"))
+ ((Vector3r,weight,Vector3r(0,0,0),,"barycentric coordinates of the projection point |yupdate|"))
((Real,relPos,0,,"position of the contact on the connection (0: node-, 1:node+) |yupdate|")),
createIndex(); /*ctor*/
);
=== added file 'pkg/common/PFacet.cpp'
--- pkg/common/PFacet.cpp 1970-01-01 00:00:00 +0000
+++ pkg/common/PFacet.cpp 2015-12-15 12:57:49 +0000
@@ -0,0 +1,648 @@
+/*****************************************************************************
+* Copyright (C) 2015 by Anna Effeindzourou anna.effeindzourou@xxxxxxxxx *
+* Copyright (C) 2015 by Bruno Chareyre bruno.chareyre@xxxxxxxxxxx *
+* Copyright (C) 2015 by Klaus Thoeni klaus.thoeni@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 "PFacet.hpp"
+#ifdef YADE_OPENGL
+ #include<lib/opengl/OpenGLWrapper.hpp>
+#endif
+
+//!################## IGeom Functors #####################
+// Function used in order to calculate the projection of the sphere on the PFacet element. The function returns if the the porjection is on the inside the triangle and the barycentric coordinates of the projection P on PFacet the element.
+boost::tuple<Vector3r,bool, double, double,double> Ig2_Sphere_PFacet_ScGridCoGeom::projection(
+ const shared_ptr<Shape>& cm2, const State& state1)
+{
+ const State* sphereSt = YADE_CAST<const State*>(&state1);
+ PFacet* Pfacet = YADE_CAST<PFacet*>(cm2.get());
+
+ vector<Vector3r> vertices;
+ vertices.push_back(Pfacet->node1->state->pos);
+ vertices.push_back(Pfacet->node2->state->pos);
+ vertices.push_back(Pfacet->node3->state->pos);
+ Vector3r center=vertices[0]+((vertices[2]-vertices[0])*(vertices[1]-vertices[0]).norm()+(vertices[1]-vertices[0])*(vertices[2]-vertices[0]).norm())/((vertices[1]-vertices[0]).norm()+(vertices[2]-vertices[1]).norm()+(vertices[0]-vertices[2]).norm());
+
+
+ Vector3r e[3] = {vertices[1]-vertices[0] ,vertices[2]-vertices[1] ,vertices[0]-vertices[2]};
+ Vector3r normal = e[0].cross(e[1])/((e[0].cross(e[1])).norm());
+
+// Vector3r centerS=sphereSt->pos+shift2;//FIXME: periodicity?
+ const Vector3r& centerS=sphereSt->pos;
+ Vector3r cl=centerS-center;
+
+ Real dist=normal.dot(cl);
+
+ if (dist<0) {normal=-normal; dist=-dist;}
+
+ Vector3r P =center+(cl - dist*normal);
+
+ Vector3r v0 = vertices[1] - vertices[0];
+ Vector3r v1 = vertices[2] - vertices[0];
+ Vector3r v2 = P - vertices[0];
+
+ // Compute dot products
+ Real dot00 = v0.dot(v0);
+ Real dot01 = v0.dot(v1);
+ Real dot02 = v0.dot(v2);
+ Real dot11 = v1.dot(v1);
+ Real dot12 = v1.dot(v2);
+
+ // Compute the barycentric coordinates of the projection P
+ Real invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
+ Real p1 = (dot11 * dot02 - dot01 * dot12) * invDenom;
+ Real p2 = (dot00 * dot12 - dot01 * dot02) * invDenom;
+ Real p3 = 1-p1-p2;
+
+ // Check if P is in triangle
+ bool isintriangle= (p1 > 0) && (p2 > 0) && (p1 + p2 < 1);
+ return boost::make_tuple(P,isintriangle,p1,p2,p3);
+}
+
+
+bool Ig2_Sphere_PFacet_ScGridCoGeom::go( const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
+ const shared_ptr<Interaction>& c)
+{
+ TIMING_DELTAS_START();
+
+ Sphere* sphere = YADE_CAST<Sphere*>(cm1.get());
+ PFacet* Pfacet = YADE_CAST<PFacet*>(cm2.get());
+
+ const State* sphereSt = YADE_CAST<const State*>(&state1);
+
+ Real sphereRadius = sphere->radius;
+ Real PFacetradius=Pfacet->radius;
+
+ vector<Vector3r> vertices;
+ vertices.push_back(Pfacet->node1->state->pos);
+ vertices.push_back(Pfacet->node2->state->pos);
+ vertices.push_back(Pfacet->node3->state->pos);
+ Vector3r center=vertices[0]+((vertices[2]-vertices[0])*(vertices[1]-vertices[0]).norm()+(vertices[1]-vertices[0])*(vertices[2]-vertices[0]).norm())/((vertices[1]-vertices[0]).norm()+(vertices[2]-vertices[1]).norm()+(vertices[0]-vertices[2]).norm());
+
+
+ Vector3r e[3] = {vertices[1]-vertices[0] ,vertices[2]-vertices[1] ,vertices[0]-vertices[2]};
+ Vector3r normal = e[0].cross(e[1])/((e[0].cross(e[1])).norm());
+
+// Vector3r centerS=sphereSt->pos+shift2;//FIXME: periodicity?
+ const Vector3r& centerS=sphereSt->pos;
+ Vector3r cl=centerS-center;
+ Real dist=normal.dot(cl);
+
+ shared_ptr<ScGridCoGeom> scm;
+ bool isNew = !(c->geom);
+
+ if (c->geom)
+ scm = YADE_PTR_CAST<ScGridCoGeom>(c->geom);
+ else
+ scm = shared_ptr<ScGridCoGeom>(new ScGridCoGeom());
+
+ if(scm->isDuplicate==2 && scm->trueInt!=c->id2)return true; //the contact will be deleted into the Law, no need to compute here.
+ scm->isDuplicate=0;
+ scm->trueInt=-1;
+
+
+ if (std::abs(dist)>(PFacetradius+sphereRadius) && !c->isReal() && !force) { // no contact, but only if there was no previous contact; ortherwise, the constitutive law is responsible for setting Interaction::isReal=false
+ TIMING_DELTAS_CHECKPOINT("Ig2_Sphere_PFacet_ScGridCoGeom");
+ return false;
+ }
+ if (dist<0) {normal=-normal; dist=-dist;}
+
+
+ boost::tuple <Vector3r,bool, double, double,double> projectionres = projection(cm2, state1);
+ Vector3r P = boost::get<0>(projectionres);
+ bool isintriangle = boost::get<1>(projectionres);
+ Real p1 = boost::get<2>(projectionres);
+ Real p2 = boost::get<3>(projectionres);
+ Real p3 = boost::get<4>(projectionres);
+
+
+ shared_ptr<Body> GridList[3]={Pfacet->conn1,Pfacet->conn2,Pfacet->conn3};
+
+ // Check if the projection of the contact point is inside the triangle
+ bool isconn1=((p1 > 0) && (p2 <= 0) && (p1 + p2 < 1))||((p1 > 0) && (p2 <= 0) && (p1 + p2 >= 1));
+ bool isconn2=((p1 > 0) && (p2 > 0) && (p1 + p2 >= 1))||((p1 <= 0) && (p2 > 0) && (p1 + p2 >= 1));
+ bool isconn3=((p1 <= 0) && (p2 > 0) && (p1 + p2 < 1))||((p1 <= 0) && (p2 <= 0) && (p1 + p2 < 1));
+
+ Real penetrationDepth=0;
+ int connnum=-1;
+
+ GridNode* GridNodeList[3]={ YADE_CAST<GridNode*>(Pfacet->node1->shape.get()), YADE_CAST<GridNode*>(Pfacet->node2->shape.get()),YADE_CAST<GridNode*>(Pfacet->node3->shape.get())};
+// If the contact projection is in the triangle, a search for an old contact is performed to export the information regarding the contact (phys)
+ if (isintriangle){
+ if(isNew){
+ for (int unsigned i=0; i<3;i++){
+ for(int unsigned j=0;j<GridNodeList[i]->pfacetList.size();j++){
+ if(GridNodeList[i]->pfacetList[j]->getId()!=c->id2){
+ boost::tuple <Vector3r,bool, double, double,double> projectionPrev = projection(GridNodeList[i]->pfacetList[j]->shape,state1);
+ bool isintrianglePrev = boost::get<1>(projectionPrev);
+ if(isintrianglePrev){ //if(!isintrianglePrev)
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,GridNodeList[i]->pfacetList[j]->getId());
+ if( intr && intr->isReal() ){// if an interaction exist between the sphere and the previous pfacet, import parameters.
+// cout<<"Copying contact with pfacet geom and phys from "<<intr->id1<<"-"<<intr->id2<<" to here ("<<c->id1<<"-"<<c->id2<<")"<<endl;
+ scm=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
+ c->geom=scm;
+ c->phys=intr->phys;
+ c->iterMadeReal=intr->iterMadeReal;// iteration from the time when the contact became real
+ scm->trueInt=c->id2;
+ scm->isDuplicate=2; //command the old contact deletion.
+ isNew=0;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ else{
+// identification of the cylinder possibly in contact with the sphere
+ if (isconn1) connnum=0;
+ if (isconn2) connnum=1;
+ if (isconn3) connnum=2;
+// check if the identified cylinder previously is in contact with the sphere
+ if (connnum!=-1){
+ //verify if there is a contact between the a neighbouring PFacet, avoid double contacts
+ for (int unsigned i=0; i<3;i++){
+ for(int unsigned j=0;j<GridNodeList[i]->pfacetList.size();j++){
+ if(GridNodeList[i]->pfacetList[j]->getId()!=c->id2){
+ boost::tuple <Vector3r,bool, double, double,double> projectionPrev = projection(GridNodeList[i]->pfacetList[j]->shape,state1);
+ bool isintrianglePrev = boost::get<1>(projectionPrev);
+ if(isintrianglePrev){
+ if (isNew){return false;}
+ else {
+// cout<<"The contact "<<c->id1<<"-"<<c->id2<<" may be copied and will be deleted now."<<endl ;
+ scm->isDuplicate=1 ;
+ scm->trueInt=-1 ;
+ return true;
+ }
+ }
+ }
+ }
+ }
+ //SPhere-cylinder contact
+ const State* sphereSt = YADE_CAST<const State*>(&state1);
+ GridConnection* gridCo = YADE_CAST<GridConnection*>(GridList[connnum]->shape.get());
+ GridNode* gridNo1 = YADE_CAST<GridNode*>(gridCo->node1->shape.get());
+ GridNode* gridNo2 = YADE_CAST<GridNode*>(gridCo->node2->shape.get());
+ State* gridNo1St = YADE_CAST<State*>(gridCo->node1->state.get());
+ State* gridNo2St = YADE_CAST<State*>(gridCo->node2->state.get());
+
+ Vector3r segt = gridCo->getSegment();
+ Real len = gridCo->getLength();
+
+ Vector3r spherePos = sphereSt->pos;
+ Vector3r branch = spherePos - gridNo1St->pos;
+ Vector3r branchN = spherePos - gridNo2St->pos;
+ for(int i=0;i<3;i++){
+ if(std::abs(branch[i])<1e-14) branch[i]=0.0;
+ if(std::abs(branchN[i])<1e-14) branchN[i]=0.0;
+ }
+ Real relPos = branch.dot(segt)/(len*len);
+ bool SGr=true;
+ if(relPos<=0){ // if the sphere projection is BEFORE the segment ...
+ if(gridNo1->ConnList.size()>1){// if the node is not an extremity of the Grid (only one connection)
+ for(int unsigned i=0;i<gridNo1->ConnList.size();i++){ // ... loop on all the Connections of the same Node ...
+ GridConnection* GC = (GridConnection*)gridNo1->ConnList[i]->shape.get();
+ if(GC==gridCo)continue;// self comparison.
+ Vector3r segtCandidate1 = GC->node1->state->pos - gridNo1St->pos; // (be sure of the direction of segtPrev to compare relPosPrev.)
+ Vector3r segtCandidate2 = GC->node2->state->pos - gridNo1St->pos;
+ Vector3r segtPrev = segtCandidate1.norm()>segtCandidate2.norm() ? segtCandidate1:segtCandidate2;
+ for(int j=0;j<3;j++){
+ if(std::abs(segtPrev[j])<1e-14) segtPrev[j]=0.0;
+ }
+ Real relPosPrev = (branch.dot(segtPrev))/(segtPrev.norm()*segtPrev.norm());
+ // ... and check whether the sphere projection is before the neighbours connections too.
+ if(relPosPrev<=0){//if the sphere projection is outside both the current Connection AND this neighbouring connection, then create the interaction if the neighbour did not already do it before.
+ for(int unsigned j=0;j<GC->pfacetList.size();j++){
+ if(GC->pfacetList[j]->getId()!=c->id2){
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,GC->pfacetList[j]->getId());
+
+ if(intr && intr->isReal()){
+ shared_ptr<ScGridCoGeom> intrGeom=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
+ if(!intrGeom->isDuplicate==1){ //skip contact.
+ if (isNew) {SGr=false;}
+ else {scm->isDuplicate=1;}/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/
+ }
+ }
+ }
+ }
+ }
+ else{//the sphere projection is outside the current Connection but inside the previous neighbour. The contact has to be handled by the Prev GridConnection, not here.
+ if (isNew) {SGr=false;}
+ else {
+ //cout<<"The contact "<<c->id1<<"-"<<c->id2<<" HAVE to be copied and deleted NOW."<<endl ;
+ scm->isDuplicate=1;
+ scm->trueInt=-1; //trueInt id de l'objet avec lequel il y a un contact -1 = rien faire
+ return true;
+ }
+ }
+ }
+ }
+ }
+
+ //Exactly the same but in the case the sphere projection is AFTER the segment.
+ else if(relPos>=1){
+ if(gridNo2->ConnList.size()>1){
+ for(int unsigned i=0;i<gridNo2->ConnList.size();i++){
+ GridConnection* GC = (GridConnection*)gridNo2->ConnList[i]->shape.get();
+ if(GC==gridCo)continue;// self comparison.
+ Vector3r segtCandidate1 = GC->node1->state->pos - gridNo2St->pos;
+ Vector3r segtCandidate2 = GC->node2->state->pos - gridNo2St->pos;
+ Vector3r segtNext = segtCandidate1.norm()>segtCandidate2.norm() ? segtCandidate1:segtCandidate2;
+ for(int j=0;j<3;j++){
+ if(std::abs(segtNext[j])<1e-14) segtNext[j]=0.0;
+ }
+
+ Real relPosNext = (branchN.dot(segtNext))/(segtNext.norm()*segtNext.norm());
+ if(relPosNext<=0){ //if the sphere projection is outside both the current Connection AND this neighbouring connection, then create the interaction if the neighbour did not already do it before.
+ for(int unsigned j=0;j<GC->pfacetList.size();j++){
+ if(GC->pfacetList[j]->getId()!=c->id2){
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,GC->pfacetList[j]->getId());
+ if(intr && intr->isReal()){
+ shared_ptr<ScGridCoGeom> intrGeom=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
+ if(!intrGeom->isDuplicate==1){
+ if (isNew) SGr=false;
+ else scm->isDuplicate=1;/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/
+ }
+ }
+ }
+ }
+ }
+ else{//the sphere projection is outside the current Connection but inside the previous neighbour. The contact has to be handled by the Prev GridConnection, not here.
+ if (isNew)SGr=false;
+ else {//cout<<"The contact "<<c->id1<<"-"<<c->id2<<" HAVE to be copied and deleted NOW."<<endl ;
+ scm->isDuplicate=1 ;
+ scm->trueInt=-1 ;
+ return true;
+ }
+ }
+ }
+ }
+ }
+ else if (relPos<=0.5){
+ if(gridNo1->ConnList.size()>1){// if the node is not an extremity of the Grid (only one connection)
+ for(int unsigned i=0;i<gridNo1->ConnList.size();i++){ // ... loop on all the Connections of the same Node ...
+ GridConnection* GC = (GridConnection*)gridNo1->ConnList[i]->shape.get();
+ if(GC==gridCo)continue;// self comparison.
+
+ Vector3r segtCandidate1 = GC->node1->state->pos - gridNo1St->pos; // (be sure of the direction of segtPrev to compare relPosPrev.)
+ Vector3r segtCandidate2 = GC->node2->state->pos - gridNo1St->pos;
+ Vector3r segtPrev = segtCandidate1.norm()>segtCandidate2.norm() ? segtCandidate1:segtCandidate2;
+ for(int j=0;j<3;j++){
+ if(std::abs(segtPrev[j])<1e-14) segtPrev[j]=0.0;
+ }
+ Real relPosPrev = (branch.dot(segtPrev))/(segtPrev.norm()*segtPrev.norm());
+ if(relPosPrev<=0){ //the sphere projection is inside the current Connection and outide this neighbour connection.
+ for(int unsigned j=0;j<GC->pfacetList.size();j++){
+ if(GC->pfacetList[j]->getId()!=c->id2){
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,GC->pfacetList[j]->getId());
+ if( intr && intr->isReal() ){// if an ineraction exist between the sphere and the previous connection, import parameters.
+ // cout<<"1Copying contact geom and phys from "<<intr->id1<<"-"<<intr->id2<<" to here ("<<c->id1<<"-"<<c->id2<<")"<<endl;
+ scm=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
+ if(isNew){
+ c->geom=scm;
+ c->phys=intr->phys;
+ c->iterMadeReal=intr->iterMadeReal;
+ }
+ scm->trueInt=c->id2;
+ scm->isDuplicate=2; //command the old contact deletion.
+ isNew=0;
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ }
+
+ else if (relPos>0.5){
+ if(gridNo2->ConnList.size()>1){
+ for(int unsigned i=0;i<gridNo2->ConnList.size();i++){
+ GridConnection* GC = (GridConnection*)gridNo2->ConnList[i]->shape.get();
+ if(GC==gridCo)continue;// self comparison.
+ Vector3r segtCandidate1 = GC->node1->state->pos - gridNo2St->pos;
+ Vector3r segtCandidate2 = GC->node2->state->pos - gridNo2St->pos;
+ Vector3r segtNext = segtCandidate1.norm()>segtCandidate2.norm() ? segtCandidate1:segtCandidate2;
+ for(int j=0;j<3;j++){
+ if(std::abs(segtNext[j])<1e-14) segtNext[j]=0.0;
+ }
+ Real relPosNext = (branchN.dot(segtNext))/(segtNext.norm()*segtNext.norm());
+ if(relPosNext<=0){ //the sphere projection is inside the current Connection and outside this neighbour connection.
+ for(int unsigned j=0;j<GC->pfacetList.size();j++){
+ if(GC->pfacetList[j]->getId()!=c->id2){
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,GC->pfacetList[j]->getId());
+ if( intr && intr->isReal() ){// if an ineraction exist between the sphere and the previous connection, import parameters.
+ scm=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
+ if(isNew){
+ c->geom=scm;
+ c->phys=intr->phys;
+ c->iterMadeReal=intr->iterMadeReal;
+ }
+ scm->trueInt=c->id2;
+ scm->isDuplicate=2; //command the old contact deletion.
+ break;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ if(SGr){
+ relPos=relPos<0?0:relPos; //min value of relPos : 0
+ relPos=relPos>1?1:relPos; //max value of relPos : 1
+ Vector3r fictiousPos=gridNo1St->pos+relPos*segt;
+ Vector3r branchF = fictiousPos - spherePos;
+
+ Real dist = branchF.norm();
+ bool SG= !(isNew && (dist > (sphere->radius + gridCo->radius)));
+ if(SG){
+ // Create the geometry :
+ if(isNew) c->geom=scm;
+ scm->radius1=sphere->radius;
+ scm->radius2=gridCo->radius;
+ scm->id3=gridCo->node1->getId();
+ scm->id4=gridCo->node2->getId();
+
+ scm->relPos=relPos;
+ Vector3r normal=branchF/dist;
+ scm->penetrationDepth = sphere->radius+gridCo->radius-dist;
+ scm->fictiousState.pos = fictiousPos;
+ scm->contactPoint = spherePos + normal*(scm->radius1 - 0.5*scm->penetrationDepth);
+ scm->fictiousState.vel = (1-relPos)*gridNo1St->vel + relPos*gridNo2St->vel;
+ scm->fictiousState.angVel =
+ ((1-relPos)*gridNo1St->angVel + relPos*gridNo2St->angVel).dot(segt/len)*segt/len //twist part : interpolated
+ + segt.cross(gridNo2St->vel - gridNo1St->vel);// non-twist part : defined from nodes velocities
+ scm->precompute(state1,scm->fictiousState,scene,c,normal,isNew,shift2,true);//use sphere-sphere precompute (with a virtual sphere)
+ return true;
+ }
+
+ }
+ }
+ else{//taking into account the shadow zone
+ for (int unsigned i=0; i<3;i++){
+ for(int unsigned j=0;j<GridNodeList[i]->pfacetList.size();j++){
+ if(GridNodeList[i]->pfacetList[j]->getId()!=c->id2){
+ boost::tuple <Vector3r,bool, double, double,double> projectionPrev = projection(GridNodeList[i]->pfacetList[j]->shape,state1);
+ bool isintrianglePrev = boost::get<1>(projectionPrev);
+ if(!isintrianglePrev){
+ //if the sphere projection is outside both the current PFacet AND this neighbouring PFacet, then create the interaction if the neighbour did not already do it before.
+ const shared_ptr<Interaction> intr = scene->interactions->find(c->id1,GridNodeList[i]->pfacetList[j]->getId());
+ if(intr && intr->isReal()){
+ shared_ptr<ScGridCoGeom> intrGeom=YADE_PTR_CAST<ScGridCoGeom>(intr->geom);
+ if(!intrGeom->isDuplicate==1){ //skip contact.
+ if (isNew) {
+ return false;}
+ else {
+ scm->isDuplicate=1;
+ }/*cout<<"Declare "<<c->id1<<"-"<<c->id2<<" as duplicated."<<endl;*/
+ }
+ }
+ }
+ else{//the sphere projection is outside the current PFacet but inside the previous neighbour. The contact has to be handled by the Prev GridConnection, not here.
+ if (isNew) {return false;}
+ else {
+ // cout<<"The contact "<<c->id1<<"-"<<c->id2<<" HAVE to be copied and deleted NOW."<<endl ;
+ scm->isDuplicate=1; //scm->isDuplicate=1;
+ scm->trueInt=-1; //trueInt id de l'objet avec lequel il y a un contact -1 = rien faire
+ return true;
+ }
+ }
+ }
+
+ }
+ }
+ }
+ }
+ if(isintriangle){
+ penetrationDepth = sphereRadius + PFacetradius - std::abs(dist);
+ normal.normalize();
+ if (penetrationDepth>0 || c->isReal() ){
+ if(isNew) c->geom=scm;
+ scm->radius1=sphereRadius;
+ scm->radius2=PFacetradius;
+ scm->id3=Pfacet->node1->getId();
+ scm->id4=Pfacet->node2->getId();
+ scm->id5=Pfacet->node3->getId();
+ scm->weight[0]=p1;
+ scm->weight[1]=p2;
+ scm->weight[2]=p3;
+ scm->penetrationDepth = penetrationDepth;
+ scm->fictiousState.pos = P;
+ scm->contactPoint = centerS - (PFacetradius-0.5*penetrationDepth)*normal;
+ scm->fictiousState.vel = (p1*Pfacet->node1->state->vel + p2*Pfacet->node2->state->vel +p3*Pfacet->node3->state->vel);
+ scm->fictiousState.angVel =
+ (p1*Pfacet->node1->state->angVel + p2*Pfacet->node2->state->angVel+ p3*Pfacet->node3->state->angVel);
+ scm->precompute(state1,scm->fictiousState,scene,c,-normal,isNew,shift2,true);//use sphere-sphere precompute (with a virtual sphere)
+ TIMING_DELTAS_CHECKPOINT("Ig2_Sphere_PFacet_ScGridCoGeom");
+ return true;
+
+ }
+ }
+ TIMING_DELTAS_CHECKPOINT("Ig2_Sphere_PFacet_ScGridCoGeom");
+ return false;
+}
+bool Ig2_Sphere_PFacet_ScGridCoGeom::goReverse( const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c)
+{
+ c->swapOrder();
+ return go(cm2,cm1,state2,state1,-shift2,force,c);
+}
+YADE_PLUGIN((Ig2_Sphere_PFacet_ScGridCoGeom));
+
+bool Ig2_GridConnection_PFacet_ScGeom::go( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c)
+{
+ GridConnection* gridCo = YADE_CAST<GridConnection*>(cm1.get());
+ PFacet* Pfacet = YADE_CAST<PFacet*>(cm2.get());
+
+ if(gridCo->node1==Pfacet->node1 || gridCo->node1==Pfacet->node2|| gridCo->node1==Pfacet->node3 || gridCo->node2==Pfacet->node1 || gridCo->node2==Pfacet->node2|| gridCo->node2==Pfacet->node3){return false;}
+
+
+
+ Body::id_t idNode1=gridCo->node1->getId();
+ Body::id_t idNode2=gridCo->node2->getId();
+ Body::id_t ids2[3]={Pfacet->conn1->getId(),Pfacet->conn2->getId(),Pfacet->conn3->getId()};
+
+ if (!scene->interactions->found(idNode1,c->id2)){
+ shared_ptr<Interaction> scm1 (new Interaction(idNode1,c->id2));
+ scene->interactions->insert(scm1);
+ }
+
+ if (!scene->interactions->found(idNode2,c->id2)){
+ shared_ptr<Interaction> scm2 (new Interaction(idNode2,c->id2));
+ scene->interactions->insert(scm2);
+ }
+
+ for (int i=0; i<3; i++){
+ int entier=i;
+ ostringstream oss;
+ string chaine = "scm";
+ oss << chaine << entier;
+ string chaine1=oss.str();
+ if (!scene->interactions->found(c->id1,ids2[i])){
+ shared_ptr<Interaction> chaine1 (new Interaction(c->id1,ids2[i]));
+ scene->interactions->insert(chaine1);
+ }
+ }
+
+ return(false);
+}
+
+bool Ig2_GridConnection_PFacet_ScGeom::goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c)
+{
+ return go(cm1,cm2,state2,state1,-shift2,force,c);
+}
+YADE_PLUGIN((Ig2_GridConnection_PFacet_ScGeom));
+
+
+bool Ig2_PFacet_PFacet_ScGeom::goOneWay( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c)
+{
+ /*FIXME : /!\ Note that this geometry doesn't take care of any unwished duplicated contact or shear force following. /!\*/
+ PFacet* Pfacet1 = YADE_CAST<PFacet*>(cm1.get());
+
+ Body::id_t idNode1=Pfacet1->node1->getId();
+ Body::id_t idNode2=Pfacet1->node2->getId();
+ Body::id_t idNode3=Pfacet1->node3->getId();
+
+ Body::id_t id2=c->id2;
+
+ if (!scene->interactions->found(idNode1,id2)){
+ shared_ptr<Interaction> scm1 (new Interaction(idNode1,id2));
+ scene->interactions->insert(scm1);
+ }
+
+ if (!scene->interactions->found(idNode2,id2)){
+ shared_ptr<Interaction> scm2 (new Interaction(idNode2,id2));
+ scene->interactions->insert(scm2);
+ }
+
+ if (!scene->interactions->found(idNode3,id2)){
+ shared_ptr<Interaction> scm3 (new Interaction(idNode3,id2));
+ scene->interactions->insert(scm3);
+ }
+
+ return (false);
+
+}
+
+
+bool Ig2_PFacet_PFacet_ScGeom::go( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c)
+{
+ PFacet* Pfacet1 = YADE_CAST<PFacet*>(cm1.get());
+ PFacet* Pfacet2 = YADE_CAST<PFacet*>(cm2.get());
+
+ if(Pfacet1->node1==Pfacet2->node1 || Pfacet1->node1==Pfacet2->node2|| Pfacet1->node1==Pfacet2->node3
+ || Pfacet1->node2==Pfacet2->node1 || Pfacet1->node2==Pfacet2->node2|| Pfacet1->node2==Pfacet2->node3
+ || Pfacet1->node3==Pfacet2->node1 || Pfacet1->node3==Pfacet2->node2|| Pfacet1->node3==Pfacet2->node3 ){ return false;}
+
+
+ goOneWay( cm1, cm2, state1, state2,shift2,force, c);
+
+ swap(c->id1,c->id2);
+
+ goOneWay( cm2, cm1, state2, state1,-shift2,force, c);
+
+ Body::id_t ids1[3]={Pfacet1->conn1->getId(),Pfacet1->conn2->getId(),Pfacet1->conn3->getId()};
+ Body::id_t ids2[3]={Pfacet2->conn1->getId(),Pfacet2->conn2->getId(),Pfacet2->conn3->getId()};
+
+ for (int i=0; i<3; i++){
+ for (int j=0; j<3; j++){
+ int entier = j+i*3;
+ ostringstream oss;
+ string chaine = "scm";
+ oss << chaine << entier;
+ string chaine1=oss.str();
+ if (!scene->interactions->found(ids1[i],ids2[j])){
+ shared_ptr<Interaction> chaine1 (new Interaction(ids1[i],ids2[j]));
+ scene->interactions->insert(chaine1);
+ }
+ }
+ }
+
+ return(false);
+}
+
+bool Ig2_PFacet_PFacet_ScGeom::goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c)
+{
+ return go(cm1,cm2,state2,state1,-shift2,force,c);
+}
+YADE_PLUGIN((Ig2_PFacet_PFacet_ScGeom));
+
+
+/********* Wall + Sphere **********/
+
+bool Ig2_Wall_PFacet_ScGeom::go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
+
+ PFacet* Pfacet = YADE_CAST<PFacet*>(cm2.get());
+
+ Body::id_t idNode1=Pfacet->node1->getId();
+ Body::id_t idNode2=Pfacet->node2->getId();
+ Body::id_t idNode3=Pfacet->node3->getId();
+
+ Body::id_t id1=c->id1;
+
+ if (!scene->interactions->found(id1,idNode1)){
+ shared_ptr<Interaction> scm1 (new Interaction(id1,idNode1));
+ scene->interactions->insert(scm1);
+ }
+
+ if (!scene->interactions->found(id1,idNode2)){
+ shared_ptr<Interaction> scm2 (new Interaction(id1,idNode2));
+ scene->interactions->insert(scm2);
+ }
+
+ if (!scene->interactions->found(id1,idNode3)){
+ shared_ptr<Interaction> scm3 (new Interaction(id1,idNode3));
+ scene->interactions->insert(scm3);
+ }
+
+ return(false);
+
+}
+
+YADE_PLUGIN((Ig2_Wall_PFacet_ScGeom));
+//!################## Bounds #####################
+
+void Bo1_PFacet_Aabb::go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r& se3, const Body* b)
+{
+ PFacet* Pfacet = YADE_CAST<PFacet*>(cm.get());
+ if(!bv){ bv=shared_ptr<Bound>(new Aabb); }
+ Aabb* aabb=static_cast<Aabb*>(bv.get());
+ Vector3r O = Pfacet->node1->state->pos;
+ Vector3r O2 =Pfacet->node2->state->pos;
+ Vector3r O3 =Pfacet->node3->state->pos;
+
+ if(!scene->isPeriodic){
+ for (int k=0;k<3;k++){
+ aabb->min[k] = min(min(O[k],O2[k]),O3[k]) - Pfacet->radius;
+ aabb->max[k] = max(max(O[k],O2[k]),O3[k]) + Pfacet->radius;
+ }
+ return;
+ }
+ else{
+ O = scene->cell->unshearPt(O);
+ O2 = scene->cell->unshearPt(O2);
+ O3= scene->cell->unshearPt(O3);
+
+ Vector3r T=scene->cell->hSize*Pfacet->cellDist.cast<Real>();
+ O = O + T;
+ O2 = O2 + T;
+ O3 = O3 + T;
+ for (int k=0;k<3;k++){
+ aabb->min[k] = min(min(O[k],O2[k]),O3[k]) - Pfacet->radius;
+ aabb->max[k] = max(max(O[k],O2[k]),O3[k]) + Pfacet->radius;
+ }
+ }
+}
+
+YADE_PLUGIN((Bo1_PFacet_Aabb));
\ No newline at end of file
=== added file 'pkg/common/PFacet.hpp'
--- pkg/common/PFacet.hpp 1970-01-01 00:00:00 +0000
+++ pkg/common/PFacet.hpp 2015-12-15 12:57:49 +0000
@@ -0,0 +1,147 @@
+/*****************************************************************************
+* Copyright (C) 2015 by Anna Effeindzourou anna.effeindzourou@xxxxxxxxx *
+* Copyright (C) 2015 by Bruno Chareyre bruno.chareyre@xxxxxxxxxxx *
+* Copyright (C) 2015 by Klaus Thoeni klaus.thoeni@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<core/Shape.hpp>
+#include <pkg/dem/ScGeom.hpp>
+#include<pkg/common/Grid.hpp>
+#include <core/Body.hpp>
+#include<pkg/common/Sphere.hpp>
+#include <pkg/common/Dispatching.hpp>
+#include <pkg/common/Grid.hpp>
+#include <pkg/dem/FrictPhys.hpp>
+#include<lib/base/Math.hpp>
+#include<pkg/common/InteractionLoop.hpp>
+#include <pkg/dem/ElasticContactLaw.hpp>
+#include <pkg/dem/Ig2_Facet_Sphere_ScGeom.hpp>
+#ifdef YADE_OPENGL
+ #include<pkg/common/GLDrawFunctors.hpp>
+#endif
+
+
+//!################## IGeom Functors ##################
+//! O/
+class Ig2_Sphere_PFacet_ScGridCoGeom: public Ig2_Sphere_GridConnection_ScGridCoGeom{
+ public :
+
+ boost::tuple<Vector3r,bool, double, double,double> projection(
+ const shared_ptr<Shape>& cm2,
+ const State& state1);
+ virtual bool go(const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c);
+ YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_PFacet_ScGridCoGeom,Ig2_Sphere_GridConnection_ScGridCoGeom,"Create/update a :yref:`ScPFaceCoGeom` instance representing intersection of :yref:`Facet` and :yref:`Sphere`.",
+ ((Real,shrinkFactor,((void)"no shrinking",0),,"The radius of the inscribed circle of the facet is decreased by the value of the sphere's radius multipled by *shrinkFactor*. From the definition of contact point on the surface made of facets, the given surface is not continuous and becomes in effect surface covered with triangular tiles, with gap between the separate tiles equal to the sphere's radius multiplied by 2×*shrinkFactor*. If zero, no shrinking is done."))
+ );
+ DECLARE_LOGGER;
+ FUNCTOR2D(Sphere,PFacet);
+ DEFINE_FUNCTOR_ORDER_2D(Sphere,PFacet);
+};
+
+REGISTER_SERIALIZABLE(Ig2_Sphere_PFacet_ScGridCoGeom);
+
+
+
+
+class Ig2_GridConnection_PFacet_ScGeom: public Ig2_Sphere_GridConnection_ScGridCoGeom{
+ public :
+ virtual bool go(const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c);
+ YADE_CLASS_BASE_DOC_ATTRS(Ig2_GridConnection_PFacet_ScGeom,Ig2_Sphere_GridConnection_ScGridCoGeom,"Create/update a :yref:`ScGeom` instance representing intersection of :yref:`Facet` and :yref:`GridConnection`.",
+ ((Real,shrinkFactor,((void)"no shrinking",0),,"The radius of the inscribed circle of the facet is decreased by the value of the sphere's radius multipled by *shrinkFactor*. From the definition of contact point on the surface made of facets, the given surface is not continuous and becomes in effect surface covered with triangular tiles, with gap between the separate tiles equal to the sphere's radius multiplied by 2×*shrinkFactor*. If zero, no shrinking is done."))
+ );
+ DECLARE_LOGGER;
+ FUNCTOR2D(GridConnection,PFacet);
+ DEFINE_FUNCTOR_ORDER_2D(GridConnection,PFacet);
+};
+
+REGISTER_SERIALIZABLE(Ig2_GridConnection_PFacet_ScGeom);
+
+//! O/
+class Ig2_PFacet_PFacet_ScGeom: public Ig2_Sphere_PFacet_ScGridCoGeom{
+ public :
+ virtual bool go(const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c);
+ virtual bool goOneWay(const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& c);
+ YADE_CLASS_BASE_DOC_ATTRS(Ig2_PFacet_PFacet_ScGeom,Ig2_Sphere_PFacet_ScGridCoGeom,"Create/update a :yref:`ScGridCoGeom` instance representing intersection of :yref:`Facet` and :yref:`Sphere`.",
+ ((Real,shrinkFactor,((void)"no shrinking",0),,"The radius of the inscribed circle of the facet is decreased by the value of the sphere's radius multipled by *shrinkFactor*. From the definition of contact point on the surface made of facets, the given surface is not continuous and becomes in effect surface covered with triangular tiles, with gap between the separate tiles equal to the sphere's radius multiplied by 2×*shrinkFactor*. If zero, no shrinking is done."))
+ );
+ DECLARE_LOGGER;
+ FUNCTOR2D(PFacet,PFacet);
+ DEFINE_FUNCTOR_ORDER_2D(PFacet,PFacet);
+};
+
+REGISTER_SERIALIZABLE(Ig2_PFacet_PFacet_ScGeom);
+
+
+/********* Wall + Sphere **********/
+
+class Ig2_Wall_PFacet_ScGeom: public Ig2_Wall_Sphere_ScGeom{
+ public:
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ YADE_CLASS_BASE_DOC_ATTRS(Ig2_Wall_PFacet_ScGeom,Ig2_Wall_Sphere_ScGeom,"Create/update a :yref:`ScGeom` instance representing intersection of :yref:`Wall` and :yref:`PFacet`.",
+ );
+ FUNCTOR2D(Wall,PFacet);
+ DEFINE_FUNCTOR_ORDER_2D(Wall,PFacet);
+};
+REGISTER_SERIALIZABLE(Ig2_Wall_PFacet_ScGeom);
+
+
+//!################## Bounds #####################
+
+class Bo1_PFacet_Aabb : public BoundFunctor
+{
+ public :
+ void go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r&, const Body*);
+ FUNCTOR1D(PFacet);
+ YADE_CLASS_BASE_DOC_ATTRS(Bo1_PFacet_Aabb,BoundFunctor,"Functor creating :yref:`Aabb` from a :yref:`PFacet`.",
+ ((Real,aabbEnlargeFactor,((void)"deactivated",-1),,"Relative enlargement of the bounding box; deactivated if negative."))
+ );
+};
+REGISTER_SERIALIZABLE(Bo1_PFacet_Aabb);
+
=== modified file 'py/gridpfacet.py'
--- py/gridpfacet.py 2015-12-11 05:43:09 +0000
+++ py/gridpfacet.py 2015-12-15 12:57:49 +0000
@@ -10,7 +10,6 @@
"""
import math,random,doctest,geom,numpy
-from yade import *
from yade.wrapper import *
try: # use psyco if available
import psyco
@@ -18,7 +17,7 @@
except ImportError: pass
from yade import utils
-from yade._utils import *
+from yade._utils import createInteraction
from minieigen import *
@@ -126,3 +125,287 @@
b.mask=mask
return b
+
+def pfacet(id1,id2,id3,wire=True,color=None,highlight=False,material=-1,mask=1,cellDist=None):
+ """
+ Create a :yref:`PFacet<PFacet>` element from 3 :yref:`GridNodes<GridNode>` which are already connected via 3 :yref:`GridConnections<GridConnection>`:
+
+ :param id1,id2,id3: already with :yref:`GridConnections<GridConnection>` connected :yref:`GridNodes<GridNode>`
+ :param bool wire: if ``True``, top and bottom facet are shown as skeleton; otherwise facets are filled.
+ :param Vector3-or-None color: color of the PFacet; random color will be assigned if ``None``.
+ :param Vector3 cellDist: for periodic boundary conditions, see :yref:`Interaction.cellDist`. Note: periodic boundary conditions are not yet implemented for PFacets!
+
+ See documentation of :yref:`yade.utils.sphere` for meaning of other parameters.
+
+ :return: Body object with the :yref:`PFacet<PFacet>` :yref:`shape<Body.shape>`.
+
+ .. note:: :yref:`GridNodes<GridNode>` and :yref:`GridConnections<GridConnection>` need to have the same radius. This is also the radius used to create the :yref:`PFacet<PFacet>`
+
+ """
+ b=Body()
+ GridN1=O.bodies[id1]; GridN2=O.bodies[id2]; GridN3=O.bodies[id3]
+ b.shape=PFacet(color=color if color else randomColor(),wire=wire,highlight=highlight,node1=GridN1,node2=GridN2,node3=GridN3)
+ GridN1.bounded=False; GridN2.bounded=False; GridN3.bounded=False
+ GridC1=O.bodies[O.interactions[id1,id2].geom.connectionBody.id]
+ GridC2=O.bodies[O.interactions[id2,id3].geom.connectionBody.id]
+ GridC3=O.bodies[O.interactions[id1,id3].geom.connectionBody.id]
+ GridC1.bounded=False
+ GridC2.bounded=False
+ GridC3.bounded=False
+
+ b.shape.conn1=GridC1
+ b.shape.conn2=GridC2
+ b.shape.conn3=GridC3
+
+ b.shape.radius=GridN1.shape.radius
+ GridC1.shape.addPFacet(b)
+ GridC2.shape.addPFacet(b)
+ GridC3.shape.addPFacet(b)
+ GridN1.shape.addPFacet(b); GridN2.shape.addPFacet(b); GridN3.shape.addPFacet(b)
+
+ V=0
+
+ utils._commonBodySetup(b,V,Vector3(0,0,0),material,pos=GridN1.state.pos,dynamic=False,fixed=True)
+ b.aspherical=False # mass and inertia are lumped into the GridNodes
+ b.mask=mask
+ return b
+
+
+#TODO: find a better way of handling the Id lists for checking duplicated gridNodes or gridNodes with the same coordinates etc. It would be better to handle this globally, maybe implement something like O.bodies.getGridNodes
+def pfacetCreator1(vertices,radius,nodesIds=[],cylIds=[],pfIds=[],wire=False,fixed=True,materialNodes=-1,material=-1,color=None):
+ """
+ Create a :yref:`PFacet<PFacet>` element from 3 vertices and automatically append to simulation. The function uses the vertices to create :yref:`GridNodes<GridNode>` and automatically checks for existing nodes.
+
+ :param [Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.
+ :param float radius: radius used to create the :yref:`PFacets<PFacet>`.
+ :param list nodesIds: list with ids of already existing :yref:`GridNodes<GridNode>`. New ids will be added.
+ :param list cylIds: list with ids of already existing :yref:`GridConnections<GridConnection>`. New ids will be added.
+ :param list pfIds: list with ids of already existing :yref:`PFacets<PFacet>`. New ids will be added.
+ :param materialNodes: specify :yref:`Body.material` of :yref:`GridNodes<GridNode>`. This material is used to make the internal connections.
+ :param material: specify :yref:`Body.material` of :yref:`PFacets<PFacet>`. This material is used for interactions with external bodies.
+
+ See documentation of :yref:`yade.utils.sphere` for meaning of other parameters.
+ """
+ n=len(nodesIds)
+ k=[0,0,0]
+ f=[0,0,0]
+ u=0
+ nod=0
+ for i in vertices:
+ u=0
+ for j in nodesIds:
+ if(i==O.bodies[j].state.pos):
+ f[nod]=j
+ k[nod]=1
+ u+=1
+ nod+=1
+ test=True
+ #if(u==0):
+ for GN in nodesIds:
+ if(i==O.bodies[GN].state.pos):
+ u=1
+ if(u==0):
+ nodesIds.append( O.bodies.append(gridNode(i,radius,wire=wire,fixed=fixed,material=materialNodes,color=color)) )
+
+ if(k==[0,0,0]):
+ pfacetCreator3(nodesIds[n],nodesIds[n+1],nodesIds[n+2],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+ if(k==[1,0,0]):
+ pfacetCreator3(f[0],nodesIds[n],nodesIds[n+1],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+ if(k==[0,1,0]):
+ pfacetCreator3(nodesIds[n],f[1],nodesIds[n+1],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+ if(k==[0,0,1]):
+ pfacetCreator3(nodesIds[n],nodesIds[n+1],f[2],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+ if(k==[1,1,0]):
+ pfacetCreator3(f[0],f[1],nodesIds[n],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+ if(k==[0,1,1]):
+ pfacetCreator3(nodesIds[n],f[1],f[2],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+ if(k==[1,0,1]):
+ pfacetCreator3(f[0],nodesIds[n],f[2],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+ if(k==[1,1,1]):
+ pfacetCreator3(f[0],f[1],f[2],cylIds=cylIds,pfIds=pfIds,wire=wire,material=material,color=color,fixed=fixed )
+
+
+def pfacetCreator2(id1,id2,vertex,radius,nodesIds=[],wire=True,materialNodes=-1,material=-1,color=None,fixed=True):
+ """
+ Create a :yref:`PFacet<PFacet>` element from 2 already existing and connected :yref:`GridNodes<GridNode>` and one vertex. The element is automatically appended to the simulation.
+
+ :param int id1,id2: ids of already with :yref:`GridConnection` connected :yref:`GridNodes<GridNode>`.
+ :param Vector3 vertex: coordinates of the vertex in the global coordinate system.
+
+ See documentation of :yref:`yade.gridpfacet.pfacetCreator1` for meaning of other parameters.
+ """
+ n=len(nodesIds)
+ nodesIds.append( O.bodies.append(gridNode(vertex,radius,wire=wire,fixed=fixed,material=materialNodes,color=color)) )
+ O.bodies.append(gridConnection(id1,nodesIds[n],radius=radius,material=materialNodes,color=color,wire=wire))
+ O.bodies.append(gridConnection(id2,nodesIds[n],radius=radius,material=materialNodes,color=color,wire=wire))
+ O.bodies.append(pfacet(id1,id2,nodesIds[n],wire=wire,material=material,color=color))
+
+
+def pfacetCreator3(id1,id2,id3,cylIds=[],pfIds=[],wire=True,material=-1,color=None,fixed=True,mask=-1):
+ """
+ Create a :yref:`PFacet` element from 3 already existing :yref:`GridNodes<GridNode>` which are not yet connected. The element is automatically appended to the simulation.
+
+ :param int id1,id2,id3: id of the 3 :yref:`GridNodes<GridNode>` forming the :yref:`PFacet`.
+
+ See documentation of :yref:`yade.gridpfacet.pfacetCreator1` for meaning of other parameters.
+ """
+ radius=O.bodies[id1].shape.radius
+ try:
+ cylIds.append(O.bodies.append(gridConnection(id1,id2,radius=radius,material=material,color=color,wire=wire,mask=mask)))
+ except:
+ pass
+ try:
+ cylIds.append(O.bodies.append(gridConnection(id2,id3,radius=radius,material=material,color=color,wire=wire,mask=mask)))
+ except:
+ pass
+ try:
+ cylIds.append(O.bodies.append(gridConnection(id3,id1,radius=radius,material=material,color=color,wire=wire,mask=mask)))
+ except:
+ pass
+ pfIds.append(O.bodies.append(pfacet(id1,id2,id3,wire=wire,material=material,color=color,mask=mask)))
+
+
+def pfacetCreator4(id1,id2,id3,pfIds=[],wire=True,material=-1,color=None,fixed=True,mask=-1):
+ """
+ Create a :yref:`PFacet<PFacet>` element from 3 already existing :yref:`GridConnections<GridConnection>`. The element is automatically appended to the simulation.
+
+ :param int id1,id2,id3: id of the 3 :yref:`GridConnections<GridConnection>` forming the :yref:`PFacet`.
+
+ See documentation of :yref:`yade.gridpfacet.pfacetCreator1` for meaning of other parameters.
+ """
+ radius=O.bodies[id1].shape.radius
+ GridN=[]
+ GridN1=O.bodies[id1].shape.node1.id
+ GridN2=O.bodies[id1].shape.node2.id
+ GridN.append(GridN1)
+ GridN.append(GridN2)
+
+ GridN1=O.bodies[id2].shape.node1.id
+ if(GridN1 not in GridN):
+ GridN.append(GridN1)
+
+ GridN2=O.bodies[id2].shape.node2.id
+ if(GridN2 not in GridN):
+ GridN.append(GridN2)
+
+ GridN1=O.bodies[id3].shape.node1.id
+ if(GridN1 not in GridN):
+ GridN.append(GridN1)
+
+ GridN2=O.bodies[id3].shape.node2.id
+ if(GridN2 not in GridN):
+ GridN.append(GridN2)
+ pfIds.append(O.bodies.append(pfacet(GridN[0],GridN[1],GridN[2],wire=wire,material=material,color=color,mask=mask)))
+
+
+def gtsPFacet(meshfile,shift=Vector3.Zero,scale=1.0,
+ radius=1,wire=True,fixed=True,materialNodes=-1,material=-1,color=None):
+ """
+ Imports mesh geometry from .gts file and automatically creates connected :yref:`PFacet3<PFacet>` elements. For an example see :ysrc:`examples/pfacet/gts-pfacet.py`.
+
+ :param string filename: .gts file to read.
+ :param [float,float,float] shift: [X,Y,Z] parameter shifts the mesh.
+ :param float scale: factor scales the mesh.
+ :param float radius: radius used to create the :yref:`PFacets<PFacet>`.
+ :param materialNodes: specify :yref:`Body.material` of :yref:`GridNodes<GridNode>`. This material is used to make the internal connections.
+ :param material: specify :yref:`Body.material` of :yref:`PFacets<PFacet>`. This material is used for interactions with external bodies.
+
+ See documentation of :yref:`yade.utils.sphere` for meaning of other parameters.
+
+ :returns: lists of :yref:`GridNode<GridNode>` ids `nodesIds`, :yref:`GridConnection<GridConnection>` ids `cylIds`, and :yref:`PFacet<PFacet>` ids `pfIds`
+ """
+ import gts,yade.pack
+ surf=gts.read(open(meshfile))
+ surf.scale(scale,scale,scale)
+ surf.translate(shift[0],shift[1],shift[2])
+ nodesIds=[]; cylIds=[]; pfIds=[]
+ for face in surf.faces():
+ a=face.vertices()[0].coords()
+ b=face.vertices()[1].coords()
+ c=face.vertices()[2].coords()
+ pfacetCreator1([a,b,c],radius=radius,nodesIds=nodesIds,cylIds=cylIds,pfIds=pfIds,wire=wire,fixed=fixed,materialNodes=materialNodes,material=material,color=color)
+ #print a,b,c
+ return nodesIds,cylIds,pfIds
+
+
+def gmshPFacet(meshfile="file.mesh",shift=Vector3.Zero,scale=1.0,orientation=Quaternion.Identity,
+ radius=1.0,wire=True,fixed=True,materialNodes=-1,material=-1,color=None):
+ """
+ Imports mesh geometry from .mesh file and automatically creates connected PFacet elements. For an example see :ysrc:`examples/pfacet/mesh-pfacet.py`.
+
+ :param string filename: .gts file to read.
+ :param [float,float,float] shift: [X,Y,Z] parameter shifts the mesh.
+ :param float scale: factor scales the mesh.
+ :param quaternion orientation: orientation of the imported geometry.
+ :param float radius: radius used to create the :yref:`PFacets<PFacet>`.
+ :param materialNodes: specify :yref:`Body.material` of :yref:`GridNodes<GridNode>`. This material is used to make the internal connections.
+ :param material: specify :yref:`Body.material` of :yref:`PFacets<PFacet>`. This material is used for interactions with external bodies.
+
+ See documentation of :yref:`yade.utils.sphere` for meaning of other parameters.
+
+ :returns: lists of :yref:`GridNode<GridNode>` ids `nodesIds`, :yref:`GridConnection<GridConnection>` ids `cylIds`, and :yref:`PFacet<PFacet>` ids `pfIds`
+
+ mesh files can easily be created with `GMSH <http://www.geuz.org/gmsh/>`_.
+
+ Additional examples of mesh-files can be downloaded from
+ http://www-roc.inria.fr/gamma/download/download.php
+ """
+ infile = open(meshfile,"r")
+ lines = infile.readlines()
+ infile.close()
+
+ nodelistVector3=[]
+ findVerticesString=0
+
+ while (lines[findVerticesString].split()[0]<>'Vertices'): # find the string with the number of Vertices
+ findVerticesString+=1
+ findVerticesString+=1
+ numNodes = int(lines[findVerticesString].split()[0])
+
+ for i in range(numNodes):
+ nodelistVector3.append(Vector3(0.0,0.0,0.0))
+ id = 0
+
+ for line in lines[findVerticesString+1:numNodes+findVerticesString+1]:
+ data = line.split()
+ nodelistVector3[id] = orientation*Vector3(float(data[0])*scale,float(data[1])*scale,float(data[2])*scale)+shift
+ id += 1
+
+ findTriangleString=findVerticesString+numNodes
+ while (lines[findTriangleString].split()[0]<>'Triangles'): # find the string with the number of Triangles
+ findTriangleString+=1
+ findTriangleString+=1
+ numTriangles = int(lines[findTriangleString].split()[0])
+
+ triList = []
+ for i in range(numTriangles):
+ triList.append([0,0,0,0])
+
+ tid = 0
+ for line in lines[findTriangleString+1:findTriangleString+numTriangles+1]:
+ data = line.split()
+ id1 = int(data[0])-1
+ id2 = int(data[1])-1
+ id3 = int(data[2])-1
+ triList[tid][0] = tid
+ triList[tid][1] = id1
+ triList[tid][2] = id2
+ triList[tid][3] = id3
+ tid += 1
+
+ nodesIds=[]; cylIds=[]; pfIds=[]
+ for i in triList:
+ a=nodelistVector3[i[1]]
+ b=nodelistVector3[i[2]]
+ c=nodelistVector3[i[3]]
+ #print 'i',i
+ #print 'a',a
+ #print 'b',b
+ #print 'c',c
+ try:
+ pfacetCreator1([a,b,c],radius=radius,nodesIds=nodesIds,cylIds=cylIds,pfIds=pfIds,
+ wire=wire,fixed=fixed,materialNodes=materialNodes,material=material,color=color)
+ except:
+ pass
+ return nodesIds,cylIds,pfIds
+
