Index: /issm/trunk/src/c/objects/Bamg/Mesh.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/Mesh.cpp (revision 5604) +++ /issm/trunk/src/c/objects/Bamg/Mesh.cpp (revision 5605) @@ -1059,53 +1059,55 @@ void Mesh::AddVertex( BamgVertex &s,Triangle* t, Icoor2* det3) { /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/Add)*/ - // ------------------------------------------- + // ------------------------------- // s2 - // ! - // /|\ ! - // / | \ ! - // / | \ ! - // tt1 / | \ tt0 ! - // / |s \ ! - // / . \ ! - // / . ` \ ! - // / . ` \ ! - // ---------------- ! + // ! + // /|\ ! + // / | \ ! + // / | \ ! + // tt1 / | \ tt0 ! + // / |s \ ! + // / . \ ! + // / . ` \ ! + // / . ` \ ! + // ---------------- ! // s0 tt2 s1 - //-------------------------------------------- - - //the three triangles tt - Triangle* tt[3]; - //three vertices of t - BamgVertex &s0 = (*t)[0], &s1=(*t)[1], &s2=(*t)[2]; - //three determinants - Icoor2 det3local[3]; - // number of zero det3 - register int nbzerodet =0; - // izerodet = egde contening the vertex s - register int izerodet=-1,iedge; + //------------------------------- + + /*Intermediaries*/ + Triangle* tt[3]; //the three triangles + Icoor2 det3local[3]; //three determinants (integer) + int nbzerodet =0; //number of zeros in det3 + int izerodet=-1; //egde containing the vertex s + int iedge; + + /*three vertices of t*/ + BamgVertex &s0=(*t)[0]; + BamgVertex &s1=(*t)[1]; + BamgVertex &s2=(*t)[2]; + //determinant of t - Icoor2 detOld = t->det; - - // infinitevertexpos = order of the infinite vertex (NULL) - // if no infinite vertex (NULL) infinitevertexpos=-1 - // else if v_i is infinite, infinitevertexpos=i - int infinitevertexpos = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0; + Icoor2 detOld=t->det; + + /* infvertexindex = index of the infinite vertex (NULL) + if no infinite vertex (NULL) infvertexindex=-1 + else if v_i is infinite, infvertexindex=i*/ + int infvertexindex = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0; //some checks - if (( infinitevertexpos <0 ) && (detOld <0) || ( infinitevertexpos >=0 ) && (detOld >0) ){ - ISSMERROR("bug in Mesh::Add, bad configuration"); + if (( infvertexindex <0 ) && (detOld <0) || ( infvertexindex >=0 ) && (detOld >0) ){ + ISSMERROR("inconsistent configuration (Contact ISSM developers)"); } // if det3 does not exist, build it - if (!det3) { + if (!det3){ //allocate det3 = det3local; //if no infinite vertex - if (infinitevertexpos<0 ) { + if (infvertexindex<0 ) { det3[0]=bamg::det(s ,s1,s2); det3[1]=bamg::det(s0,s ,s2); det3[2]=bamg::det(s0,s1,s );} else { - // one of &s1 &s2 &s0 is NULL so (&si || &sj) <=> !&sk + // one of &s1 &s2 &s0 is NULL det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ; det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ; @@ -1119,31 +1121,28 @@ //if nbzerodet>0, point s is on an egde or on a vertex - if (nbzerodet >0 ){ - if (nbzerodet == 1) { + if (nbzerodet>0){ + /*s is on an edge*/ + if (nbzerodet==1) { iedge = OppositeEdge[izerodet]; AdjacentTriangle ta = t->Adj(iedge); - // the point is on the edge - // if the point is one the boundary - // add the point in outside part - if ( t->det >=0) { // inside triangle - if ((( Triangle *) ta)->det < 0 ) { + /*if the point is one the boundary + add the point in outside part */ + if (t->det>=0){ // inside triangle + if (((Triangle*)ta)->det<0 ) { // add in outside triangle - AddVertex(s,( Triangle *) ta); + AddVertex(s,( Triangle *)ta); return; } } } - else { - //t->Echo(); - printf("\nproblem while trying to add:\n"); - s.Echo(); - ISSMERROR("Bug in Mesh::Add points duplicated %i times",nbzerodet); + else{ + ISSMERROR("Cannot add a vertex more than once. Check duplicates"); } } // remove de MarkUnSwap edge - t->SetUnMarkUnSwap(0); - t->SetUnMarkUnSwap(1); + t->SetUnMarkUnSwap(0); + t->SetUnMarkUnSwap(1); t->SetUnMarkUnSwap(2); @@ -2790,8 +2789,8 @@ * the initial simple mesh from this edge and 2 boundary triangles*/ - BamgVertex * v0=orderedvertices[0], *v1=orderedvertices[1]; + BamgVertex *v0=orderedvertices[0], *v1=orderedvertices[1]; nbt = 2; - triangles[0](0) = NULL; //infinite vertex + triangles[0](0) = NULL;//infinite vertex triangles[0](1) = v0; triangles[0](2) = v1; @@ -2823,5 +2822,4 @@ /*Now, add the vertices One by One*/ - long NbSwap=0; if (verbose>3) printf(" Begining of insertion process...\n"); @@ -2833,6 +2831,6 @@ //Find the triangle in which newvertex is located - Icoor2 dete[3]; - Triangle* tcvi = TriangleFindFromCoord(newvertex->i,dete); //(newvertex->i = integer coordinates) + Icoor2 det3[3]; + Triangle* tcvi = TriangleFindFromCoord(newvertex->i,det3); //(newvertex->i = integer coordinates) //Add newvertex to the quadtree @@ -2840,5 +2838,5 @@ //Add newvertex to the existing mesh - AddVertex(*newvertex,tcvi,dete); + AddVertex(*newvertex,tcvi,det3); //Make the mesh Delaunay around newvertex by swaping the edges @@ -2881,5 +2879,5 @@ long i; long NbSwap=0; - Icoor2 dete[3]; + Icoor2 det3[3]; //number of new points @@ -2925,5 +2923,5 @@ } vj.ReferenceNumber=0; - Triangle *tcvj=TriangleFindFromCoord(vj.i,dete); + Triangle *tcvj=TriangleFindFromCoord(vj.i,det3); if (tcvj && !tcvj->link){ tcvj->Echo(); @@ -2931,5 +2929,5 @@ } quadtree->Add(vj); - AddVertex(vj,tcvj,dete); + AddVertex(vj,tcvj,det3); NbSwap += vj.Optim(1); iv++; @@ -3569,8 +3567,8 @@ for (int icount=2; icounti,dete); + Icoor2 det3[3]; + Triangle *tcvi = TriangleFindFromCoord(vi->i,det3); quadtree->Add(*vi); - AddVertex(*vi,tcvi,dete); + AddVertex(*vi,tcvi,det3); NbSwap += vi->Optim(1,1); } @@ -4682,5 +4680,5 @@ long iv = nbvold; long NbSwap = 0; - Icoor2 dete[3]; + Icoor2 det3[3]; for (int i=nbvold;ilink) { printf("problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n"); @@ -4700,5 +4698,5 @@ ISSMERROR("!tcvi || tcvi->det < 0"); } - AddVertex(vi,tcvi,dete); + AddVertex(vi,tcvi,det3); NbSwap += vi.Optim(1); iv++; @@ -4728,5 +4726,5 @@ /*}}}1*/ /*FUNCTION Mesh::TriangleFindFromCoord{{{1*/ -Triangle * Mesh::TriangleFindFromCoord(const I2 & B,Icoor2 dete[3], Triangle *tstart) const { +Triangle * Mesh::TriangleFindFromCoord(const I2 & B,Icoor2 det3[3], Triangle *tstart) const { /*Original code from Frederic Hecht (BAMG v1.01, Mesh2.cpp/FindTriangleContening)*/ @@ -4769,7 +4767,7 @@ ISSMASSERT(k0>=0);// k0 the NULL vertex int k1=NextVertex[k0],k2=PreviousVertex[k0]; - dete[k0]=det(B,(*t)[k1],(*t)[k2]); - dete[k1]=dete[k2]=-1; - if (dete[k0] > 0) // outside B + det3[k0]=det(B,(*t)[k1],(*t)[k2]); + det3[k1]=det3[k2]=-1; + if (det3[k0] > 0) // outside B return t; t = t->TriangleAdj(OppositeEdge[k0]); @@ -4787,15 +4785,15 @@ j= OppositeVertex[jj]; - dete[j] = detop; //det(*b,*s1,*s2); + det3[j] = detop; //det(*b,*s1,*s2); jn = NextVertex[j]; jp = PreviousVertex[j]; - dete[jp]= det(*(*t)(j),*(*t)(jn),B); - dete[jn] = t->det-dete[j] -dete[jp]; - - // count the number k of dete <0 + det3[jp]= det(*(*t)(j),*(*t)(jn),B); + det3[jn] = t->det-det3[j] -det3[jp]; + + // count the number k of det3 <0 int k=0,ii[3]; - if (dete[0] < 0 ) ii[k++]=0; - if (dete[1] < 0 ) ii[k++]=1; - if (dete[2] < 0 ) ii[k++]=2; + if (det3[0] < 0 ) ii[k++]=0; + if (det3[1] < 0 ) ii[k++]=1; + if (det3[2] < 0 ) ii[k++]=2; // 0 => ok // 1 => go in way 1 @@ -4803,8 +4801,6 @@ if (k==0) break; - if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]); - if ( k>=3){ - ISSMERROR("k>=3"); - } + if (k==2 && BinaryRand()) Exchange(ii[0],ii[1]); + ISSMASSERT(k<3); AdjacentTriangle t1 = t->Adj(jj=ii[0]); if ((t1.det() < 0 ) && (k == 2)) @@ -4812,10 +4808,10 @@ t=t1; j=t1;// for optimisation we now the -det[OppositeVertex[j]]; - detop = -dete[OppositeVertex[jj]]; + detop = -det3[OppositeVertex[jj]]; jj = j; } if (t->det<0) // outside triangle - dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop; + det3[0]=det3[1]=det3[2]=-1,det3[OppositeVertex[jj]]=detop; return t; } Index: /issm/trunk/src/c/objects/Bamg/Triangle.cpp =================================================================== --- /issm/trunk/src/c/objects/Bamg/Triangle.cpp (revision 5604) +++ /issm/trunk/src/c/objects/Bamg/Triangle.cpp (revision 5605) @@ -259,7 +259,11 @@ /*FUNCTION Triangle::SetAdj2{{{1*/ void Triangle::SetAdj2(short a,Triangle *t,short aat){ - adj[a]=t; //the adjacent triangle to the edge a is t - AdjEdgeIndex[a]=aat; //position of the edge in the adjacent triangle - if(t) { //if t!=NULL add adjacent triangle to t (this) + /*For current triangle: + * - a is the index of the edge were the adjency is set (in [0 2]) + * - t is the adjacent triangle + * - aat is the index of the same edge in the adjacent triangle*/ + adj[a]=t; + AdjEdgeIndex[a]=aat; + if(t){ //if t!=NULL add adjacent triangle to t (this) t->adj[aat]=this; t->AdjEdgeIndex[aat]=a;