Index: /issm/trunk/src/c/Bamgx/meshtype.h =================================================================== --- /issm/trunk/src/c/Bamgx/meshtype.h (revision 3233) +++ /issm/trunk/src/c/Bamgx/meshtype.h (revision 3234) @@ -59,7 +59,4 @@ template inline T Min3 (const T &a,const T & b,const T & c){return Min(Min(a,b),c);} - //some functions - Int4 AGoodNumberPrimeWith(Int4 n); - //Inline functions inline int BinaryRand(){ Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp (revision 3233) +++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp (revision 3234) @@ -5,4 +5,5 @@ #include "../BamgObjects.h" +#include "../shared/shared.h" #undef __FUNCT__ Index: /issm/trunk/src/c/Bamgx/objects/Triangle.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangle.cpp (revision 3233) +++ /issm/trunk/src/c/Bamgx/objects/Triangle.cpp (revision 3234) @@ -5,4 +5,5 @@ #include "../BamgObjects.h" +#include "../shared/shared.h" #undef __FUNCT__ Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3233) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3234) @@ -5,4 +5,5 @@ #include "../BamgObjects.h" +#include "../shared/shared.h" #include "QuadTree.h" #include "SetOfE4.h" Index: /issm/trunk/src/c/Bamgx/objects/Vertex.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 3233) +++ /issm/trunk/src/c/Bamgx/objects/Vertex.cpp (revision 3234) @@ -5,4 +5,5 @@ #include "../BamgObjects.h" +#include "../shared/shared.h" #undef __FUNCT__ Index: /issm/trunk/src/c/Bamgx/shared/AGoodNumberPrimeWith.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/shared/AGoodNumberPrimeWith.cpp (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/AGoodNumberPrimeWith.cpp (revision 3234) @@ -0,0 +1,29 @@ +#include "AGoodNumberPrimeWith.h" + +Int4 AGoodNumberPrimeWith(Int4 n){ + + //list of big prime numbers + const Int4 BigPrimeNumber[] ={ 567890359L, + 567890431L, 567890437L, 567890461L, 567890471L, + 567890483L, 567890489L, 567890497L, 567890507L, + 567890591L, 567890599L, 567890621L, 567890629L , 0}; + + //initialize o and pi + Int4 o =0; + Int4 pi=BigPrimeNumber[1]; + + //loop until BigPrimeNumber[i]==0 (end of BigPrimeNumber) + for (int i=0; BigPrimeNumber[i]; i++){ + + //compute r, rest of the remainder of the division of BigPrimeNumber[i] by n + Int4 r = BigPrimeNumber[i] % n; + + /*compute oo = min ( r , n-r , |n - 2r|, |n-3r|)*/ + Int4 oo = Min(Min(r,n-r),Min(Abs(n-2*r),Abs(n-3*r))); + if ( o < oo){ + o=oo; + pi=BigPrimeNumber[i]; + } + } + return pi; +} Index: /issm/trunk/src/c/Bamgx/shared/AGoodNumberPrimeWith.h =================================================================== --- /issm/trunk/src/c/Bamgx/shared/AGoodNumberPrimeWith.h (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/AGoodNumberPrimeWith.h (revision 3234) @@ -0,0 +1,11 @@ +#ifndef _FINDTRIANGLE_H_ +#define _FINDTRIANGLE_H_ + +#include "../meshtype.h" + +namespace bamg { + + Int4 AGoodNumberPrimeWith(Int4 n); + +} +#endif Index: /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdge.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdge.cpp (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdge.cpp (revision 3234) @@ -0,0 +1,38 @@ +#include "CloseBoundaryEdge.h" + +TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) { + int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; + int dir=0; + if (k<0){ + throw ErrorException(__FUNCT__,exprintf("k<0")); + } + int kkk=0; + Icoor2 IJ_IA,IJ_AJ; + TriangleAdjacent edge(t,OppositeEdge[k]); + for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) { + kkk++; + if (kkk>=1000){ + throw ErrorException(__FUNCT__,exprintf("kkk>=1000")); + } + Vertex &vI = *edge.EdgeVertex(0); + Vertex &vJ = *edge.EdgeVertex(1); + I2 I=vI, J=vJ, IJ= J-I; + IJ_IA = (IJ ,(A-I)); + if (IJ_IA<0) { + if (dir>0) {a=1;b=0;return edge;}// change of signe => I + else {dir=-1; + continue;}};// go in direction i + IJ_AJ = (IJ ,(J-A)); + if (IJ_AJ<0) { + if(dir<0) {a=0;b=1;return edge;} + else {dir = 1; + continue;}}// go in direction j + double IJ2 = IJ_IA + IJ_AJ; + if (IJ2==0){ + throw ErrorException(__FUNCT__,exprintf("IJ2==0")); + } + a= IJ_AJ/IJ2; + b= IJ_IA/IJ2; + return edge; + } +} Index: /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdge.h =================================================================== --- /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdge.h (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdge.h (revision 3234) @@ -0,0 +1,16 @@ +#ifndef _CLOSEBOUNDARYEDGE_H_ +#define _CLOSEBOUNDARYEDGE_H_ + +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../toolkits/toolkits.h" + +#include "../objects/Triangle.h" + +namespace bamg { + + TriangleAdjacent CloseBoundaryEdge(I2 ,Triangle *, double &,double &) ; + +} +#endif Index: /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdgeV2.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdgeV2.cpp (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/CloseBoundaryEdgeV2.cpp (revision 3234) @@ -0,0 +1,148 @@ +#include "CloseBoundaryEdgeV2.h" + +TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) { + // walk around the vertex + // version 2 for remove the probleme if we fill the hole + //int bug=1; + // Triangle *torigine = t; + // restart: + // int dir=0; + if (t->link != 0){ + throw ErrorException(__FUNCT__,exprintf("t->link != 0")); + } + // to have a starting edges + // try the 3 edge bourna-- in case of internal hole + // and choice the best + // + // the probleme is in case of the fine and long internal hole + // for exemple neart the training edge of a wing + Vertex * s=0,*s1=0, *s0=0; + Icoor2 imax = MaxICoor22; + Icoor2 l0 = imax,l1 = imax; + double dd2 = imax;// infinity + TriangleAdjacent er; + int cas=-2; + for (int j=0;j<3;j++) + { + TriangleAdjacent ta=t->FindBoundaryEdge(j); + if (! (Triangle *) ta) continue; + s0 = ta.EdgeVertex(0); + s1 = ta.EdgeVertex(1); + I2 A = * s0; + I2 B = *ta.EdgeVertex(1); + I2 AB = B-A,AC=C-A,BC=B-C; + Icoor2 ACAC = (AC,AC), BCBC = (BC,BC); + Icoor2 AB2 = Norme2_2(AB); // ||AB||^2 + Icoor2 ABAC = (AB,AC); // AB.AC| + + double d2; + if ( ABAC < 0 ) // DIST A + { + if ( (d2=(double) ACAC) < dd2) + { + er = ta; + l0 = ACAC; + l1 = BCBC; + cas = 0; + s = s0; + } + } + else if (ABAC > AB2) // DIST B + { + if ( (d2=(double) BCBC) < dd2) + { + dd2 = d2; + er = Adj(ta); // other direction + l0 = BCBC; + l1 = ACAC; + cas = 1; + s = s1; + } + } + else // DIST AB + { + + double det_2 = (double) Det(AB,AC); + det_2 *= det_2; // square of area*2 of triangle ABC + d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC + + if (d2 < dd2) + { + dd2 = d2; + er = ta; + l0 = (AC,AC); + l1 = (BC,BC); + s = 0; + cas = -1; + b = ((double) ABAC/(double) AB2); + a = 1 - b; + } + } + } + if (cas ==-2){ + throw ErrorException(__FUNCT__,exprintf("cas==-2")); + } + // l1 = ||C s1|| , l0 = ||C s0|| + // where s0,s1 are the vertex of the edge er + + if ( s) + { + t=er; + TriangleAdjacent edge(er); + + int kkk=0; + int linkp = t->link == 0; + + Triangle * tt=t=edge=Adj(Previous(edge)); + do { // loop over vertex s + kkk++; + if (edge.EdgeVertex(0)!=s && kkk>=10000){ + throw ErrorException(__FUNCT__,exprintf("edge.EdgeVertex(0)!=s && kkk>=10000")); + } + + int link = tt->link == 0; + if ((link + linkp) == 1) + { // a boundary edge + Vertex * st = edge.EdgeVertex(1); + I2 I=*st; + Icoor2 ll = Norme2_2 (C-I); + if (ll < l1) { // the other vertex is neart + s1=st; + l1=ll; + er = edge; + if(lllink) + { ta = ta.Adj(); + tc=ta; + k = ta; + Exchange(aa,bb); + if (!tc->link){ + throw ErrorException(__FUNCT__,exprintf("!tc->link")); + } + } + a[VerticesOfTriangularEdge[k][0]] = aa; + a[VerticesOfTriangularEdge[k][1]] = bb; + a[OppositeVertex[k]] = 1- aa -bb; + return Th.Number(tc); + } +} Index: /issm/trunk/src/c/Bamgx/shared/FindTriangle.h =================================================================== --- /issm/trunk/src/c/Bamgx/shared/FindTriangle.h (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/FindTriangle.h (revision 3234) @@ -0,0 +1,19 @@ +#ifndef _FINDTRIANGLE_H_ +#define _FINDTRIANGLE_H_ + +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../toolkits/toolkits.h" + +#include "../meshtype.h" +#include "../objects/Triangle.h" + +namespace bamg { + +// class Triangles; + + Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside); + +} +#endif Index: /issm/trunk/src/c/Bamgx/shared/ForceEdge.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/shared/ForceEdge.cpp (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/ForceEdge.cpp (revision 3234) @@ -0,0 +1,74 @@ +#include "ForceEdge.h" + +int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) { + int NbSwap =0; + if (!a.t || !b.t){ // the 2 vertex is in a mesh + throw ErrorException(__FUNCT__,exprintf("!a.t || !b.t")); + } + int k=0; + taret=TriangleAdjacent(0,0); // erreur + + TriangleAdjacent tta(a.t,EdgesVertexTriangle[a.vint][0]); + Vertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2; + // we turn around a in the direct sens + + Icoor2 det2 = v2 ? det(*v2,a,b): -1 , det1; + if(v2) // normal case + det2 = det(*v2,a,b); + else { // no chance infini vertex try the next + tta= Previous(Adj(tta)); + v2 = tta.EdgeVertex(0); + vbegin =v2; + if (!v2){ + throw ErrorException(__FUNCT__,exprintf("!v2")); + } + det2 = det(*v2,a,b); + } + + while (v2 != &b) { + TriangleAdjacent tc = Previous(Adj(tta)); + v1 = v2; + v2 = tc.EdgeVertex(0); + det1 = det2; + det2 = v2 ? det(*v2,a,b): det2; + + if((det1 < 0) && (det2 >0)) { + // try to force the edge + Vertex * va = &a, *vb = &b; + tc = Previous(tc); + if (!v1 || !v2){ + throw ErrorException(__FUNCT__,exprintf("!v1 || !v2")); + } + Icoor2 detss = 0,l=0,ks; + while ((ks=SwapForForcingEdge( va, vb, tc, detss, det1,det2,NbSwap))) + if(l++ > 10000000) { + throw ErrorException(__FUNCT__,exprintf("Loop in forcing Egde, nb de swap=%i, nb of try swap (%i) too big",NbSwap,l)); + } + Vertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1); + if (( aa == &a ) && (bb == &b) || (bb == &a ) && (aa == &b)) { + tc.SetLock(); + a.Optim(1,0); + b.Optim(1,0); + taret = tc; + return NbSwap; + } + else + { + taret = tc; + return -2; // error boundary is crossing + } + } + tta = tc; + k++; + if (k>=2000){ + throw ErrorException(__FUNCT__,exprintf("k>=2000")); + } + if ( vbegin == v2 ) return -1;// error + } + + tta.SetLock(); + taret=tta; + a.Optim(1,0); + b.Optim(1,0); + return NbSwap; +} Index: /issm/trunk/src/c/Bamgx/shared/ForceEdge.h =================================================================== --- /issm/trunk/src/c/Bamgx/shared/ForceEdge.h (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/ForceEdge.h (revision 3234) @@ -0,0 +1,20 @@ +#ifndef _FORCEEDGE_H_ +#define _FORCEEDGE_H_ + +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../toolkits/toolkits.h" + +#include "../objects/Vertex.h" +#include "../objects/TriangleAdjacent.h" + +namespace bamg { + +// class Vertex; +// class TriangleAdjacent; + + int ForceEdge(Vertex &a, Vertex & b,TriangleAdjacent & taret) ; + +} +#endif Index: /issm/trunk/src/c/Bamgx/shared/SwapForForcingEdge.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/shared/SwapForForcingEdge.cpp (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/SwapForForcingEdge.cpp (revision 3234) @@ -0,0 +1,87 @@ +#include "SwapForForcingEdge.h" + +int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb ,TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & NbSwap) { + // l'arete ta coupe l'arete pva pvb + // de cas apres le swap sa coupe toujours + // on cherche l'arete suivante + // on suppose que detsa >0 et detsb <0 + // attention la routine echange pva et pvb + + if(tt1.Locked()) return 0; // frontiere croise + + TriangleAdjacent tt2 = Adj(tt1); + Triangle *t1=tt1,*t2=tt2;// les 2 triangles adjacent + Int1 a1=tt1,a2=tt2;// les 2 numero de l arete dans les 2 triangles + if ( a1<0 || a1>=3 ){ + throw ErrorException(__FUNCT__,exprintf("a1<0 || a1>=3")); + } + + Vertex & sa= (* t1)[VerticesOfTriangularEdge[a1][0]]; + Vertex & s1= (*t1)[OppositeVertex[a1]]; + Vertex & s2= (*t2)[OppositeVertex[a2]]; + + + Icoor2 dets2 = det(*pva,*pvb,s2); + Icoor2 det1=t1->det , det2=t2->det ; + Icoor2 detT = det1+det2; + if ((det1<=0 ) || (det2<=0)){ + throw ErrorException(__FUNCT__,exprintf("(det1<=0 ) || (det2<=0)")); + } + if ( (detsa>=0) || (detsb<=0) ){ // [a,b] cut infinite line va,bb + throw ErrorException(__FUNCT__,exprintf("(detsa>=0) || (detsb<=0)")); + } + Icoor2 ndet1 = bamg::det(s1,sa,s2); + Icoor2 ndet2 = detT - ndet1; + + int ToSwap =0; //pas de swap + if ((ndet1 >0) && (ndet2 >0)) + { // on peut swaper + if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0)) + ToSwap =1; + else // swap alleatoire + if (BinaryRand()) + ToSwap =2; + } + if (ToSwap) NbSwap++, + bamg::swap(t1,a1,t2,a2,&s1,&s2,ndet1,ndet2); + + int ret=1; + + if (dets2 < 0) {// haut + dets1 = ToSwap ? dets1 : detsa ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if (dets2 > 0){// bas + dets1 = ToSwap ? dets1 : detsb ; + detsb = dets2; + //xxxx tt1 = ToSwap ? tt1 : Next(tt2); + if(!ToSwap) tt1 = Next(tt2); + } + else { // changement de sens + ret = -1; + Exchange(pva,pvb); + Exchange(detsa,detsb); + Exchange(dets1,dets2); + Exchange(tt1,tt2); + dets1=-dets1; + dets2=-dets2; + detsa=-detsa; + detsb=-detsb; + + if (ToSwap) + if (dets2 < 0) {// haut + dets1 = (ToSwap ? dets1 : detsa) ; + detsa = dets2; + tt1 = Previous(tt2) ;} + else if (dets2 > 0){// bas + dets1 = (ToSwap ? dets1 : detsb) ; + detsb = dets2; + if(!ToSwap) tt1 = Next(tt2); + } + else {// on a fin ??? + tt1 = Next(tt2); + ret =0;} + + } + return ret; +} Index: /issm/trunk/src/c/Bamgx/shared/SwapForForcingEdge.h =================================================================== --- /issm/trunk/src/c/Bamgx/shared/SwapForForcingEdge.h (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/SwapForForcingEdge.h (revision 3234) @@ -0,0 +1,21 @@ +#ifndef _SWAPFORFORCINGEDGE_H_ +#define _SWAPFORFORCINGEDGE_H_ + +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../toolkits/toolkits.h" + +#include "../meshtype.h" +#include "../objects/Vertex.h" +#include "../objects/TriangleAdjacent.h" + +namespace bamg { + +// class Vertex; +// class TriangleAdjacent; + + int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , TriangleAdjacent & tt1,Icoor2 & dets1, Icoor2 & detsa,Icoor2 & detsb, int & nbswap); + +} +#endif Index: /issm/trunk/src/c/Bamgx/shared/shared.h =================================================================== --- /issm/trunk/src/c/Bamgx/shared/shared.h (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/shared.h (revision 3234) @@ -0,0 +1,18 @@ +/* \file shared.h + * \brief: header file for all shared routines. + */ + + +#ifndef _SHAREDBamg_H_ +#define _SHAREDBamg_H_ + +#include "AGoodNumberPrimeWith.h" +#include "CloseBoundaryEdge.h" +#include "CloseBoundaryEdgeV2.h" +#include "FindTriangle.h" +#include "ForceEdge.h" +#include "shared.h" +#include "SwapForForcingEdge.h" +#include "swap.h" + +#endif Index: /issm/trunk/src/c/Bamgx/shared/swap.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/shared/swap.cpp (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/swap.cpp (revision 3234) @@ -0,0 +1,44 @@ +#include "swap.h" + +void swap(Triangle *t1,Int1 a1, Triangle *t2,Int1 a2, Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2){ + // -------------------------------------------------------------- + // Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles + // + // sb sb + // / | \ / \ ! + // as1/ | \ /a2 \ ! + // / | \ / t2 \ ! + // s1 /t1 | t2 \s2 --> s1 /___as2___\s2 ! + // \ a1|a2 / \ as1 / + // \ | / \ t1 / + // \ | / as2 \ a1/ + // \ | / \ / + // sa sa + // ------------------------------------------------------------- + int as1 = NextEdge[a1]; + int as2 = NextEdge[a2]; + int ap1 = PreviousEdge[a1]; + int ap2 = PreviousEdge[a2]; + (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb + (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa + // mise a jour des 2 adjacences externes + TriangleAdjacent taas1 = t1->Adj(as1), + taas2 = t2->Adj(as2), + tas1(t1,as1), tas2(t2,as2), + ta1(t1,a1),ta2(t2,a2); + // externe haut gauche + taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap()); + // externe bas droite + taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap()); + // remove the Mark UnMarkSwap + t1->SetUnMarkUnSwap(ap1); + t2->SetUnMarkUnSwap(ap2); + // interne + tas1.SetAdj2(tas2); + + t1->det = det1; + t2->det = det2; + + t1->SetTriangleContainingTheVertex(); + t2->SetTriangleContainingTheVertex(); +} Index: /issm/trunk/src/c/Bamgx/shared/swap.h =================================================================== --- /issm/trunk/src/c/Bamgx/shared/swap.h (revision 3234) +++ /issm/trunk/src/c/Bamgx/shared/swap.h (revision 3234) @@ -0,0 +1,21 @@ +#ifndef _SWAP_H_ +#define _SWAP_H_ + +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../toolkits/toolkits.h" + +#include "../meshtype.h" +#include "../objects/Vertex.h" +#include "../objects/Triangle.h" + +namespace bamg { + + //class Vertex; + //class Triangle; + + void swap(Triangle *t1,Int1 a1, Triangle *t2,Int1 a2, Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2); + +} +#endif