yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06770
[Branch ~yade-dev/yade/trunk] Rev 2670: - update periodicity documentation
------------------------------------------------------------
revno: 2670
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Fri 2011-01-21 20:55:58 +0100
message:
- update periodicity documentation
- rename Hsize->hSize to be consistent with names convention.
- remove a duplicated ref and add a CGAL ref
modified:
core/Cell.cpp
core/Cell.hpp
doc/references.bib
doc/sphinx/formulation.rst
pkg/common/InteractionLoop.cpp
pkg/dem/HertzMindlin.cpp
pkg/dem/PeriIsoCompressor.cpp
pkg/dem/Shop.cpp
pkg/dem/VTKRecorder.cpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'core/Cell.cpp'
--- core/Cell.cpp 2011-01-20 14:55:15 +0000
+++ core/Cell.cpp 2011-01-21 19:55:58 +0000
@@ -7,13 +7,13 @@
// total transformation; M = (Id+G).M = F.M
trsf+=_trsfInc*trsf;
invTrsf=trsf.inverse();
- // Hsize contains colums with updated base vectors
- Hsize+=_trsfInc*Hsize;
- if(Hsize.determinant()==0){ throw runtime_error("Cell is degenerate (zero volume)."); }
+ // hSize contains colums with updated base vectors
+ hSize+=_trsfInc*hSize;
+ if(hSize.determinant()==0){ throw runtime_error("Cell is degenerate (zero volume)."); }
// lengths of transformed cell vectors, skew cosines
Matrix3r Hnorm; // normalized transformed base vectors
for(int i=0; i<3; i++){
- Vector3r base(Hsize.col(i));
+ Vector3r base(hSize.col(i));
_size[i]=base.norm(); base/=_size[i]; //base is normalized now
Hnorm(0,i)=base[0]; Hnorm(1,i)=base[1]; Hnorm(2,i)=base[2];};
// skew cosines
@@ -27,7 +27,7 @@
// pure unshear transformation
_unshearTrsf=_shearTrsf.inverse();
// some parts can branch depending on presence/absence of shear
- _hasShear=(Hsize(0,1)!=0 || Hsize(0,2)!=0 || Hsize(1,0)!=0 || Hsize(1,2)!=0 || Hsize(2,0)!=0 || Hsize(2,1)!=0);
+ _hasShear=(hSize(0,1)!=0 || hSize(0,2)!=0 || hSize(1,0)!=0 || hSize(1,2)!=0 || hSize(2,0)!=0 || hSize(2,1)!=0);
// OpenGL shear matrix (used frequently)
fillGlShearTrsfMatrix(_glShearTrsfMatrix);
=== modified file 'core/Cell.hpp'
--- core/Cell.hpp 2011-01-20 14:55:15 +0000
+++ core/Cell.hpp 2011-01-21 19:55:58 +0000
@@ -82,23 +82,23 @@
Real norm=x/sz; period=(int)floor(norm); return (norm-period)*sz;}
// relative position and velocity for interaction accross multiple cells
- Vector3r intrShiftPos(const Vector3i& cellDist) const { return Hsize*cellDist.cast<Real>(); }
- Vector3r intrShiftVel(const Vector3i& cellDist) const { if(homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND) return velGrad*Hsize*cellDist.cast<Real>(); return Vector3r::Zero(); }
+ Vector3r intrShiftPos(const Vector3i& cellDist) const { return hSize*cellDist.cast<Real>(); }
+ Vector3r intrShiftVel(const Vector3i& cellDist) const { if(homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND) return velGrad*hSize*cellDist.cast<Real>(); return Vector3r::Zero(); }
// return body velocity while taking away mean field velocity (coming from velGrad) if the mean field velocity is applied on velocity
Vector3r bodyFluctuationVel(const Vector3r& pos, const Vector3r& vel) const { if(homoDeform==HOMO_VEL || homoDeform==HOMO_VEL_2ND) return (vel-velGrad*pos); return vel; }
Vector3r getRefSize() const { return refSize; }
- // change length of initial cell vectors (scale current Hsize by the ratio)
- // should be identical to: { refSize=s; Hsize=invTrsf*Hsize; for(int i=0;i<3;i++) Hsize.col(i)=Hsize.col(i).normalized()*refSize[i]; /* transform back*/ Hsize=trsf*Hsize; integrateAndUpdate(0); }
- void setRefSize(const Vector3r& s){ for(int i=0;i<3;i++) Hsize.col(i)*=s[i]/refSize[i]; refSize=s; integrateAndUpdate(0); }
- Matrix3r getHsize() const { return Hsize; }
- void setHsize(const Matrix3r& m){ Hsize=m; for(int k=0;k<3;k++) refSize[k]=Hsize.col(k).norm(); integrateAndUpdate(0);}
- Matrix3r getHsize0() const { return invTrsf*Hsize; }
+ // change length of initial cell vectors (scale current hSize by the ratio)
+ // should be identical to: { refSize=s; hSize=invTrsf*hSize; for(int i=0;i<3;i++) hSize.col(i)=hSize.col(i).normalized()*refSize[i]; /* transform back*/ hSize=trsf*hSize; integrateAndUpdate(0); }
+ void setRefSize(const Vector3r& s){ for(int i=0;i<3;i++) hSize.col(i)*=s[i]/refSize[i]; refSize=s; integrateAndUpdate(0); }
+ Matrix3r getHSize() const { return hSize; }
+ void setHSize(const Matrix3r& m){ hSize=m; for(int k=0;k<3;k++) refSize[k]=hSize.col(k).norm(); integrateAndUpdate(0);}
+ Matrix3r getHSize0() const { return invTrsf*hSize; }
Matrix3r getTrsf() const { return trsf; }
- void setTrsf(const Matrix3r& m){ Hsize=m*invTrsf*Hsize; trsf=m; integrateAndUpdate(0); }
+ void setTrsf(const Matrix3r& m){ hSize=m*invTrsf*hSize; trsf=m; integrateAndUpdate(0); }
// return current cell volume
- Real getVolume() const {return Hsize.determinant();}
+ Real getVolume() const {return hSize.determinant();}
void postLoad(Cell&){ integrateAndUpdate(0); }
// to resolve overloads
@@ -106,25 +106,26 @@
Vector3r wrapPt_py(const Vector3r& pt) const { return wrapPt(pt);}
enum { HOMO_NONE=0, HOMO_POS=1, HOMO_VEL=2, HOMO_VEL_2ND=3 };
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Cell,Serializable,"Parameters of periodic boundary conditions. Only applies if O.isPeriodic==True.",
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(Cell,Serializable,"Parameters of periodic boundary conditions. Only applies if O.isPeriodic==True.",
/* overridden below to be modified by getters/setters because of intended side-effects */
((Vector3r,refSize,Vector3r(1,1,1),,"[overridden]"))
((Matrix3r,trsf,Matrix3r::Identity(),,"[overridden]"))
- ((Matrix3r,Hsize,Matrix3r::Identity(),,"[overridden]"))
+ ((Matrix3r,hSize,Matrix3r::Identity(),,"[overridden]"))
/* normal attributes */
((Matrix3r,invTrsf,Matrix3r::Identity(),Attr::readonly,"Inverse of current transformation matrix of the cell."))
((Matrix3r,velGrad,Matrix3r::Zero(),,"Velocity gradient of the transformation; used in :yref:`NewtonIntegrator`."))
((Matrix3r,prevVelGrad,Matrix3r::Zero(),Attr::readonly,"Velocity gradient in the previous step."))
((int,homoDeform,3,Attr::triggerPostLoad,"Deform (:yref:`velGrad<Cell.velGrad>`) the cell homothetically, by adjusting positions or velocities of particles. The values have the following meaning: 0: no homothetic deformation, 1: set absolute particle positions directly (when ``velGrad`` is non-zero), but without changing their velocity, 2: adjust particle velocity (only when ``velGrad`` changed) with Δv_i=Δ ∇v x_i. 3: as 2, but include a 2nd order term in addition -- the derivative of 1 (convective term in the velocity update).")),
+ /*deprec*/ ((Hsize,hSize,"conform to usual DEM terminology")),
+ /*init*/,
/*ctor*/ integrateAndUpdate(0),
/*py*/
// override some attributes above
- .add_property("refSize",&Cell::getRefSize,&Cell::setRefSize,"Reference size of the cell (lengths of initial cell vectors, i.e. columns of :yref:`Hsize<Cell.Hsize>`). Modifying this value will scale cell vectors such to have desired length.")
- .add_property("Hsize",&Cell::getHsize,&Cell::setHsize,"Base cell vectors (columns of the matrix), updated at every step from :yref:`velGrad<Cell.velGrad>` (:yref:`trsf<Cell.trsf>` accumulates applied :yref:`velGrad<Cell.velGrad>` transformations).")
+ .add_property("refSize",&Cell::getRefSize,&Cell::setRefSize,"Reference size of the cell (lengths of initial cell vectors, i.e. columns of :yref:`hSize<Cell.hSize>`). Modifying this value will scale cell vectors such to have desired length.")
+ .add_property("hSize",&Cell::getHSize,&Cell::setHSize,"Base cell vectors (columns of the matrix), updated at every step from :yref:`velGrad<Cell.velGrad>` (:yref:`trsf<Cell.trsf>` accumulates applied :yref:`velGrad<Cell.velGrad>` transformations).")
.add_property("trsf",&Cell::getTrsf,&Cell::setTrsf,"Current transformation matrix of the cell.")
// useful properties
- .add_property("Hsize0",&Cell::getHsize0,"Initial value of Hsize, before applying transformation (computed as :yref:`invTrsf<Cell.invTrsf>` × :yref:`Hsize<Cell.Hsize>`.")
+ .add_property("hSize0",&Cell::getHSize0,"Initial value of hSize, before applying transformation (computed as :yref:`invTrsf<Cell.invTrsf>` × :yref:`hSize<Cell.hSize>`.")
.def_readonly("size",&Cell::getSize_copy,"Current size of the cell, i.e. lengths of 3 cell lateral vectors after applying current trsf. Update automatically at every step.")
.add_property("volume",&Cell::getVolume,"Current volume of the cell.")
// debugging only
=== modified file 'doc/references.bib'
--- doc/references.bib 2010-12-09 05:46:51 +0000
+++ doc/references.bib 2011-01-21 19:55:58 +0000
@@ -384,10 +384,22 @@
doi = "10.1680/geot.1979.29.1.47"
}
-@Article{ Hassan2010,
- author = "A. Hassan and B. Chareyre and F. Darve and J. Meyssonier and F. Flin",
- title = "Microtomography-based Discrete Element Modelling of Creep in Snow",
- journal = "Granular Matter",
- year = "2010 (submitted)"
+@Article{ CundallStrack1979,
+ author = "P.A. Cundall and O.D.L. Strack",
+ title = "A discrete numerical model for granular assemblies",
+ journal = "Geotechnique",
+ volume = "",
+ number = "29",
+ pages = "47--65",
+ year = "1979",
+ doi = "10.1680/geot.1979.29.1.47"
}
+@Article{ cgal,
+ Author = "Jean-Daniel Boissonnat and Olivier Devillers and Sylvain Pion and Monique Teillaud and Mariette Yvinec",
+ Title = "Triangulations in CGAL",
+ Journal = "Computational Geometry: Theory and Applications",
+ Volume = "22",
+ Pages = "5--19",
+ Year = "2002"
+}
=== modified file 'doc/sphinx/formulation.rst'
--- doc/sphinx/formulation.rst 2011-01-21 08:14:28 +0000
+++ doc/sphinx/formulation.rst 2011-01-21 19:55:58 +0000
@@ -723,7 +723,7 @@
\vec{K}_{ij}=\sum_k (K_{Nk}-K_{Tk})\vec{n}_{i}\vec{n}_{j}+K_{Tk}=\sum_j K_{Nk}\left((1-\xi)\vec{n}_{i}\vec{n}_{j}+\xi\right)
-with $i$ and $j\in\{x,y,z\}$. Equations :eq:`eq-dtcr-axes` and :eq:`eq-dtcr-particle-stiffness` determine $\Dtcr$ in a simulation. A similar approach generalized to all 6 DOFs is implemented by the :yref:`GlobalStiffnessTimeStepper` engine in Yade. The derivation of generalized stiffness including rotational terms is very similar but not developed here, for simplicity. For full reference, see "PFC3D - Theoretical Background".
+with $i$ and $j\in\{x,y,z\}$. Equations :eq:`eq-dtcr-axes` and :eq:`eq-dtcr-particle-stiffness` determine $\Dtcr$ in a simulation. A similar approach generalized to all 6 DOFs is implemented by the :yref:`GlobalStiffnessTimeStepper` engine in Yade. The derivation of generalized stiffness including rotational terms is very similar but not developped here, for simplicity. For full reference, see "PFC3D - Theoretical Background".
Note that for computation efficiency reasons, eigenvalues of the stiffness matrices are not computed. They are only approximated assuming than DOF's are uncoupled, and using diagonal terms of $K.M^{-1}$. They give good approximates in typical mechanical systems.
@@ -835,18 +835,27 @@
Periodic boundary conditions
============================
-While most DEM simulations happen in $R^3$ space, it is frequently useful to avoid boundary effects by using periodic space instead. In order to satisfy periodicity conditions, periodic space is created by repetition of parallelepiped-shaped cell. In Yade, periodic space is implemented in the :yref:`Cell` class. The cell is determined by
-
-* \item the reference size $\vec{s}$ (:yref:`Cell.refSize`), giving reference cell configuration (which is always perpendicular): axis-aligned cuboid with corners $(0,0,0)$ and $\vec{s}$;
-* \item the transformation matrix $\mat{T}$ (:yref:`Cell.trsf`).
-
-The transformation matrix $\mat{T}$ can hold arbitrary linear transformation composed of scaling, rotation and shear. Volume change of the cell is given by $\det\mat{T}$. The cell can be manipulated by directly changing its transformation matrix $\mat{T}$ and its reference size $\vec{s}$.
-
-Additionally, we define teransformation gradient $\nabla \vec{v}$ (:yref:`Cell.velGrad`) which can be automatically integrated at every step using the Euler scheme
-
-.. math:: \next{\mat{T}}=\curr{\mat{T}}+\nabla \vec{v} \Dt.
-
-Along with the automatic integration of cell transformation, there is an option to homothetically displace all particles so that $\nabla \vec{v}$ is swept linearly over the whole simulation (enabled via :yref:`NewtonIntegrator.homotheticCellResize`). This avoids all boundary effects coming from change of the transformation.
+While most DEM simulations happen in $R^3$ space, it is frequently useful to avoid boundary effects by using periodic space instead. In order to satisfy periodicity conditions, periodic space is created by repetition of parallelepiped-shaped cell. In Yade, periodic space is implemented in the :yref:`Cell` class. The geometry of the cell in the reference coordinates system is defined by three edges of the parallepiped. The corresponding base vectors are stored in the columns of matrix :yref:`Cell.hSize`.
+
+The initial :yref:`Cell.hSize` can be explicitely defined as a 3x3 matrix at the begining of the simulation. There are no restricitions on the possible shapes: any parallelepiped is accepted as the initial undeformed period.
+
+The deformation of the period over time is defined via a matrix representing the gradient of an homogeneous velocity field $\nabla \vec{v}$ (:yref:`Cell.velGrad`). This gradient represents arbitrary combinations of rotations and stretches. It can be imposed externaly or updated by periodic engines in order to reach target strain values or to maintain some prescribed stress.
+The velocity gradient is integrated automatically over time, and the cumulated transformation is reflected in the transformation matrix $\mat{F}$ (:yref:`Cell.trsf`). :yref:`Cell.hSize` will also be updated. The update reads (it is similar for hSize), with $I$ the identity matrix:
+
+.. math:: \next{\mat{F}}=(I+\nabla \vec{v} \Dt)\curr{\mat{F}}.
+
+There is an alternative way to define the initial period geometry, using:
+
+* \ the reference size $\vec{s}$ (:yref:`Cell.refSize`), defining the dimensions or a rectangular parallelepiped with corners $(0,0,0)$ and $\vec{s}$;
+* \ the transformation matrix $\mat{F}$ (:yref:`Cell.trsf`).
+
+For instance, if :yref:`Cell.refSize` is set equal to the unit vector and if :yref:`Cell.trsf` is defined as a rotation matrix, the initial period will be an inclined cube. This method and the previous one (direct assignement of hSize) are not exactly equivalent: assigning hSize keeps the transformation :yref:`Cell.trsf` unaffected, hence the simulation will start with a null deformation (:yref:`Cell.trsf` =I); inversely, if the period geometry is defined using :yref:`Cell.trsf`$=m$, then the initial value will be $m$.
+
+It is believed that the first method is generally more convenient, since it will let :yref:`Cell.trsf` reflect only the transformation produced during the simulation, independently of the initial period geometry.
+
+In all cases, the period geometry should not be modified during a simulation, be it via hSize, trsf, or refSize. The velocity gradient :yref:`Cell.velGrad` is the only variable that let the period deformation be correctly accounted for in constitutive laws and Newton integrator (:yref:`NewtonIntegrator`).
+
+Along with the automatic integration of cell transformation, there is an option to homothetically displace all particles so that $\nabla \vec{v}$ is applied over the whole simulation (enabled via :yref:`Cell.homoDeform`). This avoids all boundary effects coming from change of the transformation.
Collision detection in periodic cell
------------------------------------
=== modified file 'pkg/common/InteractionLoop.cpp'
--- pkg/common/InteractionLoop.cpp 2010-12-13 12:11:43 +0000
+++ pkg/common/InteractionLoop.cpp 2011-01-21 19:55:58 +0000
@@ -53,7 +53,7 @@
size_t callbacksSize=callbacks.size();
// cache transformed cell size
- Matrix3r cellHsize; if(scene->isPeriodic) cellHsize=scene->cell->Hsize;
+ Matrix3r cellHsize; if(scene->isPeriodic) cellHsize=scene->cell->hSize;
// force removal of interactions that were not encountered by the collider
// (only for some kinds of colliders; see comment for InteractionContainer::iterColliderLastRun)
=== modified file 'pkg/dem/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp 2011-01-17 18:21:10 +0000
+++ pkg/dem/HertzMindlin.cpp 2011-01-21 19:55:58 +0000
@@ -196,8 +196,8 @@
Vector3r shearIncrement;
if(nonLin>1){
State *de1=Body::byId(id1,scene)->state.get(), *de2=Body::byId(id2,scene)->state.get();
- Vector3r shiftVel=scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->Hsize*contact->cellDist.cast<Real>()) : Vector3r::Zero();
- Vector3r shift2 = scene->isPeriodic ? Vector3r(scene->cell->Hsize*contact->cellDist.cast<Real>()): Vector3r::Zero();
+ Vector3r shiftVel=scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->hSize*contact->cellDist.cast<Real>()) : Vector3r::Zero();
+ Vector3r shift2 = scene->isPeriodic ? Vector3r(scene->cell->hSize*contact->cellDist.cast<Real>()): Vector3r::Zero();
Vector3r incidentV = geom->getIncidentVel(de1, de2, scene->dt, shift2, shiftVel, /*preventGranularRatcheting*/ nonLin>2 );
=== modified file 'pkg/dem/PeriIsoCompressor.cpp'
--- pkg/dem/PeriIsoCompressor.cpp 2011-01-19 18:22:56 +0000
+++ pkg/dem/PeriIsoCompressor.cpp 2011-01-21 19:55:58 +0000
@@ -99,7 +99,7 @@
for (int i=0;i<3;i++) strain[i]=log(scene->cell->trsf(i,i));
//Compute volume of the deformed cell
- Real volume=scene->cell->Hsize.determinant();
+ Real volume=scene->cell->hSize.determinant();
//Compute sum(fi*lj) and stiffness
stressTensor = Matrix3r::Zero();
@@ -114,7 +114,7 @@
GenericSpheresContact* gsc=YADE_CAST<GenericSpheresContact*> ( I->geom.get() );
//Contact force
Vector3r f= ( reversedForces?-1.:1. ) * ( nsi->normalForce+nsi->shearForce );
- Vector3r branch=Body::byId(I->getId2(),scene)->state->pos + scene->cell->Hsize*I->cellDist.cast<Real>() -Body::byId(I->getId1(),scene)->state->pos;
+ Vector3r branch=Body::byId(I->getId2(),scene)->state->pos + scene->cell->hSize*I->cellDist.cast<Real>() -Body::byId(I->getId1(),scene)->state->pos;
stressTensor+=f*branch.transpose();
if( !dynCell ){
for ( int i=0; i<3; i++ ) sumStiff[i]+=abs ( gsc->normal[i] ) *nsi->kn+ ( 1-abs ( gsc->normal[i] ) ) *nsi->ks;
@@ -218,7 +218,7 @@
}
for (int k=0;k<3;k++) strainRate[k]=scene->cell->velGrad(k,k);
//Update energy input
- Real dW=(scene->cell->velGrad*stressTensor).trace()*scene->dt*scene->cell->Hsize.determinant();
+ Real dW=(scene->cell->velGrad*stressTensor).trace()*scene->dt*scene->cell->hSize.determinant();
externalWork+=dW;
if(scene->trackEnergy) scene->energy->add(-dW,"velGradWork",velGradWorkIx,/*non-incremental*/false);
prevGrow = strainRate;
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2011-01-19 18:22:56 +0000
+++ pkg/dem/Shop.cpp 2011-01-21 19:55:58 +0000
@@ -577,7 +577,7 @@
void Shop::getStressLWForEachBody(vector<Matrix3r>& bStresses, bool revertSign){
const shared_ptr<Scene>& scene=Omega::instance().getScene();
bStresses.resize(scene->bodies->size());
- for (int k=0;k<scene->bodies->size();k++) bStresses[k]=Matrix3r::Zero();
+ for (unsigned int k=0;k<scene->bodies->size();k++) bStresses[k]=Matrix3r::Zero();
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
GenericSpheresContact* geom=YADE_CAST<GenericSpheresContact*>(I->geom.get());
@@ -589,7 +589,7 @@
if (!scene->isPeriodic)
bStresses[I->getId2()]+=(3.0/(4.0*Mathr::PI*pow(geom->refR2,3)))*f*((geom->contactPoint- (Body::byId(I->getId2(),scene)->state->pos)).transpose());
else
- bStresses[I->getId2()]+=(3.0/(4.0*Mathr::PI*pow(geom->refR2,3)))*f* ((geom->contactPoint- (Body::byId(I->getId2(),scene)->state->pos + (scene->cell->Hsize*I->cellDist.cast<Real>()))).transpose());
+ bStresses[I->getId2()]+=(3.0/(4.0*Mathr::PI*pow(geom->refR2,3)))*f* ((geom->contactPoint- (Body::byId(I->getId2(),scene)->state->pos + (scene->cell->hSize*I->cellDist.cast<Real>()))).transpose());
}
}
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2010-12-21 12:19:50 +0000
+++ pkg/dem/VTKRecorder.cpp 2011-01-21 19:55:58 +0000
@@ -230,8 +230,8 @@
// create two line objects; each of them has one endpoint inside the cell and the other one sticks outside
// A,B are the "fake" bodies outside the cell for id1 and id2 respectively, p1,p2 are the displayed points
// distance in cell units for shifting A away from p1; negated value is shift of B away from p2
- Vector3r ptA(p01+scene->cell->Hsize*(wrapCellDist[I->getId2()]-I->cellDist).cast<Real>());
- Vector3r ptB(p02+scene->cell->Hsize*(wrapCellDist[I->getId1()]-I->cellDist).cast<Real>());
+ Vector3r ptA(p01+scene->cell->hSize*(wrapCellDist[I->getId2()]-I->cellDist).cast<Real>());
+ Vector3r ptB(p02+scene->cell->hSize*(wrapCellDist[I->getId1()]-I->cellDist).cast<Real>());
vtkIdType idPtA=intrBodyPos->InsertNextPoint(ptA[0],ptA[1],ptA[2]), idPtB=intrBodyPos->InsertNextPoint(ptB[0],ptB[1],ptB[2]);
vtkSmartPointer<vtkLine> line1B(vtkSmartPointer<vtkLine>::New()); line1B->GetPointIds()->SetId(0,I->getId1()); line1B->GetPointIds()->SetId(1,idPtB);
vtkSmartPointer<vtkLine> lineA2(vtkSmartPointer<vtkLine>::New()); lineA2->GetPointIds()->SetId(0,idPtA); line1B->GetPointIds()->SetId(1,I->getId2());
