Commit 7354a84f authored by Thomas Pitiot 's avatar Thomas Pitiot
Browse files

working on particules

parent 981277ad
......@@ -39,14 +39,18 @@ public :
typedef typename PFP::MAP MAP;
typedef typename PFP::VEC3 VEC3;
typedef VertexAttribute<VEC3, MAP> TAB_POS;
typedef FaceAttribute<VEC3,MAP> TAB_FACE;
typedef VolumeAttribute<VEC3,MAP> TAB_VOL;
MAP& m;
const TAB_POS& position;
const TAB_FACE& face_center;
const TAB_VOL& volume_center;
VEC3 m_positionFace;
VEC3 m_positionVolume;
unsigned int crossCell ;
......@@ -54,13 +58,17 @@ public :
ParticleCell3D(MAP& map) : m(map)
{}
ParticleCell3D(MAP& map, Dart belonging_cell, VEC3 pos, const TAB_POS& tabPos) :
ParticleCell3D(MAP& map, Dart belonging_cell, VEC3 pos, const TAB_POS& tabPos,const TAB_FACE& fa_center = NULL,
const TAB_VOL& vol_center=NULL) :
Algo::MovingObjects::ParticleBase<PFP>(pos),
m(map),
position(tabPos),
d(belonging_cell)
{
m_positionFace = pointInFace(d);
if(vol_center) volume_center(vol_center);
if(fa_center) face_center(fa_center);
reset_positionFace();
reset_positionVolume();
this->setState(VOLUME);
}
......@@ -76,23 +84,13 @@ public :
d = cell;
}
inline VEC3 pointInFace(Dart d);
inline Geom::Orientation3D isLeftENextVertex(const VEC3& c, Dart d, const VEC3& base);
inline bool isRightVertex(const VEC3& c, Dart d, const VEC3& base);
inline Geom::Orientation3D whichSideOfFace(const VEC3& c, Dart d);
inline Geom::Orientation3D isLeftL1DVol(const VEC3& c, Dart d, const VEC3& base, const VEC3& top);
inline Geom::Orientation3D isRightDVol(const VEC3& c, Dart d, const VEC3& base, const VEC3& top);
inline Geom::Orientation3D isAbove(const VEC3& c, Dart d, const VEC3& top);
// inline int isLeftL1DFace(const VEC3& c, Dart d, const VEC3& base, const VEC3& normal);
// inline int isRightDFace(const VEC3& c, Dart d, const VEC3& base, const VEC3& normal);
inline Geom::Orientation3D whichSideOfPlanVolume(const VEC3& c, Dart d, const VEC3& base, const VEC3& top);
inline int whichSideOfPlan(const VEC3& c, Dart d, const VEC3& base, const VEC3& normal); // orientation par rapport au plan de gauche de l'arete visée
......@@ -114,9 +112,10 @@ public :
void volumeState(const VEC3& current);
void volumeSpecialCase(const VEC3& current);
void placeOnRightFaceAndRightEdge(const VEC3& current,bool * casON, bool * enDessous);
void reset_positionFace();
void reset_positionVolume();
void move(const VEC3& newCurrent)
{
......
......@@ -52,8 +52,6 @@ public :
void volumeState(const VEC3& current, CellMarkerMemo<MAP, VOLUME>& memo_cross);
void volumeSpecialCase(const VEC3& current, CellMarkerMemo<MAP, VOLUME>& memo_cross);
std::vector<Dart> move(const VEC3& newCurrent, CellMarkerMemo<MAP, VOLUME>& memo_cross);
std::vector<Dart> move(const VEC3& newCurrent);
......
......@@ -139,6 +139,7 @@ void ParticleCell3DMemo<PFP>::vertexState(const VEC3& current, CellMarkerMemo<MA
} while(verif && this->d!=ddd);
if(verif) {
this->reset_positionVolume();
volumeState(current,memo_cross);
return;
}
......@@ -179,8 +180,9 @@ void ParticleCell3DMemo<PFP>::vertexState(const VEC3& current, CellMarkerMemo<MA
std::cout << "numerical rounding ?" << std::endl;
this->d = dd;
this->m_position = this->pointInFace(this->d);
this->m_positionFace = this->m_position;
this->reset_positionFace();
this->reset_positionVolume();
this->m_position = this->m_positionFace ;
volumeState(current,memo_cross);
return;
}
......@@ -190,9 +192,10 @@ void ParticleCell3DMemo<PFP>::vertexState(const VEC3& current, CellMarkerMemo<MA
if(wsof != 0)
{
this->m_position = this->pointInFace(this->d);
this->d = this->nextNonPlanar(this->d);
this->m_positionFace = this->pointInFace(this->d);
this->reset_positionFace();
this->reset_positionVolume();
this->m_position = this->m_positionFace ;
volumeState(current,memo_cross);
}
else
......@@ -232,6 +235,7 @@ void ParticleCell3DMemo<PFP>::edgeState(const VEC3& current, CellMarkerMemo<MAP,
switch (this->whichSideOfEdge(current,this->d)) {
case Geom::OVER : // on est sur la face en question
this->d=this->m.phi1(this->d);
this->reset_positionFace();
faceState(current,memo_cross);
return;
case Geom::ON :// on est sur l'arete, il faut tester les sommets
......@@ -273,8 +277,9 @@ void ParticleCell3DMemo<PFP>::edgeState(const VEC3& current, CellMarkerMemo<MAP,
// d = nextNonPlanar(m.phi1(d));
this->d=this->m.phi1(this->d);
this->m_position = this->pointInFace(this->d);
this->m_positionFace = this->m_position;
this->reset_positionFace();
this->reset_positionVolume();
this->m_position = this->m_positionFace ;
volumeState(current,memo_cross);
return;
}
......@@ -298,7 +303,8 @@ void ParticleCell3DMemo<PFP>::edgeState(const VEC3& current, CellMarkerMemo<MAP,
#endif
this->d = this->m.phi3(this->d);
this->d = this->nextNonPlanar(this->d);
this->m_positionFace = this->pointInFace(this->d);
this->reset_positionFace();
this->reset_positionVolume();
volumeState(current, memo_cross);
return;
}
......@@ -307,7 +313,8 @@ void ParticleCell3DMemo<PFP>::edgeState(const VEC3& current, CellMarkerMemo<MAP,
std::cout << "wsof UNDER" << std::endl;
#endif
this->d = this->nextNonPlanar(this->d);
this->m_positionFace = this->pointInFace(this->d);
this->reset_positionFace();
this->reset_positionVolume();
volumeState(current, memo_cross);
return;
}
......@@ -400,236 +407,58 @@ void ParticleCell3DMemo<PFP>::volumeState(const VEC3& current, CellMarkerMemo<MA
std::cout << "volumeStateMemo " << this->d << " " << this->m_position<< " " << this->m_positionFace<< std::endl;
#endif
if(!memo_cross.isMarked(this->d)) memo_cross.mark(this->d);
CellMarkerStore<MAP, FACE> mark(this->m);
bool above;
Geom::Orientation3D testRight=Geom::OVER;
do {
above=false;
Dart dd=this->d;
bool particularcase=false;
Geom::Orientation3D testLeft = this->isLeftL1DVol(current,this->d,this->m_positionFace,this->m_position);
if(testLeft!=Geom::UNDER) {
this->d = this->m.phi_1(this->d);
while(dd!=this->d && this->isLeftL1DVol(current,this->d,this->m_positionFace,this->m_position)!=Geom::UNDER)
this->d = this->m.phi_1(this->d);
if(dd==this->d)
particularcase=true;
}
else {
testRight = this->isRightDVol(current,this->d,this->m_positionFace,this->m_position);
while(testRight!=Geom::OVER && dd!=this->m.phi1(this->d)) {
this->d = this->m.phi1(this->d);
testRight = this->isRightDVol(current,this->d,this->m_positionFace,this->m_position);
}
if(testLeft==0 && dd==this->m.phi1(this->d))
particularcase=true;
}
if(particularcase) //(this->m_position,this->m_positionFace,c) presque alignés et si c est proche de la face
//aucun des "above" sur les dart ne va donner de résultats concluant (proche de 0 pour tous)
{
if(isnan(current[0]) || isnan(current[1]) || isnan(current[2]))
{
std::cout << __FILE__ << " " << __LINE__ << " NaN !" << std::endl;
this->m_position = current;
return;
}
#ifdef DEBUG
std::cout << "particularCase" << std::endl;
#endif
this->volumeSpecialCase(current, memo_cross);
return;
}
Geom::Orientation3D testAbove = this->isAbove(current,this->d,this->m_position);
if(testAbove!=Geom::UNDER || (testRight==Geom::ON && this->isAbove(current,this->m.phi_1(this->d),this->m_position)!=Geom::UNDER)) {
if(testAbove==Geom::OVER || this->whichSideOfFace(current,this->d)==Geom::UNDER) {
mark.mark(this->d);
above=true;
this->d = this->m.phi2(this->d);
if(mark.isMarked(this->d)) {
dd = this->d;
this->d = this->nextFaceNotMarked(this->d,mark);
mark.mark(this->d);
if(this->d==dd) {
#ifdef DEBUG
std::cout << "dd=d" << std::endl;
#endif
this->volumeSpecialCase(current,memo_cross);
return;
}
}
this->m_positionFace = this->pointInFace(this->d);
}
}
} while(above);
Geom::Orientation3D wsof = this->whichSideOfFace(current,this->d);
bool casON=false;
bool enDessous=false;
this->placeOnRightFaceAndRightEdge(current,&casON,&enDessous);
Geom::Orientation3D wsof =this->whichSideOfFace(current,this->d);
switch (wsof) {
case Geom::OVER :
this->m_position = current;
this->setState(VOLUME);
break;
case Geom::ON : // on est sur la face
this->m_position = this->m_positionFace = current;
if(enDessous)// on teste par rapport a la derniere face du tetra si en dessous on est sur la face, sinon c'est qu'on était "on" et on est sur l'edge
{
this->setState(FACE);
}
else
{
if(casON)
{
this->setState(VERTEX);
}
else
{
this->setState(EDGE);
}
if(wsof==Geom::UNDER) {
this->m_position = current;
this->setState(VOLUME);
}
else if(wsof==Geom::ON) {
if(this->isAbove(current,this->d,this->m_position)==Geom::UNDER) {
this->m_position = this->m_positionFace = current;
this->setState(FACE);
}
else {
this->m_position = this->m_positionFace = current;
edgeState(current, memo_cross);
}
}
else {
if(this->isAbove(current,this->d,this->m_position)==Geom::UNDER) {
// this->m_position = m.intersectDartPlaneLine(d,this->m_position,current);
Algo::Surface::Geometry::intersectionLineConvexFace<PFP>(this->m,this->d,this->position,this->m_position,current-this->m_position,this->m_position);
faceState(current, memo_cross,wsof);
}
else {
// this->m_position = m.intersectDartPlaneLineEdge(d,this->m_position,current);
}
return;
default : // under l'obj est de l'autre coté de la face
Algo::Surface::Geometry::intersectionLineConvexFace<PFP>(this->m,this->d,this->position,this->m_position,current-this->m_position,this->m_position);
edgeState(current, memo_cross);
}
}
if(enDessous)
{
this->reset_positionFace();
faceState(current,memo_cross,wsof);
}
else
{
if(casON)
{
vertexState(current,memo_cross);
}
else
{
edgeState(current,memo_cross);
}
}
return;
}
}
template <typename PFP>
void ParticleCell3DMemo<PFP>::volumeSpecialCase(const VEC3& current, CellMarkerMemo<MAP, VOLUME>& memo_cross)
{
#ifdef DEBUG
std::cout << "volumeSpecialCaseMemo " << this->d << " " << this->m_position<< " " << this->m_positionFace<< std::endl;
#endif
if(!memo_cross.isMarked(this->d)) memo_cross.mark(this->d);
Dart dd;
CellMarkerStore<MAP, FACE> mark(this->m);
Dart d_min;
std::vector<Dart> dart_list;
std::vector<float> dist_list;
std::list<Dart> visitedFaces; // Faces that are traversed
visitedFaces.push_back(this->d); // Start with the face of d
std::list<Dart>::iterator face;
// For every face added to the list
for (face = visitedFaces.begin();face != visitedFaces.end(); ++face)
{
if (!mark.isMarked(*face))
{ // Face has not been visited yet
Geom::Orientation3D wsof = this->whichSideOfFace(current,*face);
if(wsof==Geom::OVER)
{
this->d = *face;
if(this->isAbove(current,this->d,this->m_position)==Geom::UNDER)
{
// this->m_position = m.intersectDartPlaneLine(d,this->m_position,current);
Algo::Surface::Geometry::intersectionLineConvexFace<PFP>(this->m,this->d,this->position,this->m_position,current-this->m_position,this->m_position);
faceState(current, memo_cross,Geom::OVER);
}
else
{
// this->m_position = m.intersectDartPlaneLineEdge(d,this->m_position,current);
Algo::Surface::Geometry::intersectionLineConvexFace<PFP>(this->m,this->d,this->position,this->m_position,current-this->m_position,this->m_position);
edgeState(current, memo_cross);
}
return;
}
else if(wsof==Geom::ON)
{
this->m_position = current;
this->d = *face;
faceState(current, memo_cross);
return;
}
Geom::Plane3D<typename PFP::REAL> pl = Algo::Surface::Geometry::facePlane<PFP>(this->m,*face,this->position);
if(pl.normal()*VEC3(current-this->m_position)>0)
{
dist_list.push_back(-pl.distance(current));
dart_list.push_back(*face);
}
// If the face wasn't crossed then mark visited darts (current face)
// and add non visited adjacent faces to the list of face
Dart ddd = *face;
do {
mark.mark(ddd);
Dart adj = this->m.phi2(ddd); // Get adjacent face
if (adj != ddd && !mark.isMarked(adj))
visitedFaces.push_back(adj);// Add it
} while(ddd!=*face);
}
}
if(dist_list.size()>0) {
float min=dist_list[0];
for(unsigned int i = 1;i<dist_list.size();++i)
{
if(dist_list[i]<min)
{
this->d=dart_list[i];
min = dist_list[i];
}
}
}
this->m_positionFace = this->pointInFace(this->d);
Geom::Orientation3D wsof = this->whichSideOfFace(current,this->d);
if(wsof==Geom::UNDER) {
this->m_position = current;
this->setState(VOLUME);
}
else if(wsof==Geom::ON) {
if(this->isAbove(current,this->d,this->m_position)==Geom::UNDER) {
this->m_position = current;
this->setState(FACE);
}
else {
this->m_position = current;
this->setState(EDGE);
}
}
else {
if(this->isAbove(current,this->d,this->m_position)==Geom::UNDER)
{
// this->m_position = m.intersectDartPlaneLine(d,this->m_position,current);
Algo::Surface::Geometry::intersectionLineConvexFace<PFP>(this->m,this->d,this->position,this->m_position,current-this->m_position,this->m_position);
faceState(current, memo_cross,wsof);
}
else {
// this->m_position = m.intersectDartPlaneLineEdge(d,this->m_position,current);
Algo::Surface::Geometry::intersectionLineConvexFace<PFP>(this->m,this->d,this->position,this->m_position,current-this->m_position,this->m_position);
edgeState(current, memo_cross);
}
}
}
} // namespace MovingObjects
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment