Index: /issm/trunk/src/c/Bamgx/Bamgx.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2876) +++ /issm/trunk/src/c/Bamgx/Bamgx.cpp (revision 2877) @@ -84,6 +84,5 @@ //build metric if not given in input if (verbosity>1) printf(" Generating Metric...\n"); - for(i=0;i no continuite GeometricalEdge * Adj [2]; - int SensAdj[2]; + int DirAdj[2]; // private: int flag ; @@ -656,5 +656,4 @@ int static counter; // to kown the number of mesh in memory - int OnDisk; // true if on disk Geometry & Gh; // Geometry Triangles & BTh; // Background Mesh Bth==*this =>no background @@ -713,15 +712,16 @@ Triangles(BamgMesh* bamgmesh,BamgOpts* bamgopts); - Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) - :Gh(BT.Gh),BTh(BT) { - try {GeomToTriangles1(nbvx,keepBackVertices);} - catch(...) { this->~Triangles(); throw; } } - - Triangles(Int4 nbvx,Geometry & G) - :Gh(G),BTh(*this){ - try { GeomToTriangles0(nbvx);} - catch(...) { this->~Triangles(); throw; } } - Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR - Triangles(const Triangles &,const int *flag,const int *bb); // truncature + Triangles(Int4 nbvx,Triangles & BT,int keepBackVertices=1) :Gh(BT.Gh),BTh(BT) { + try {GeomToTriangles1(nbvx,keepBackVertices);} + catch(...) { this->~Triangles(); throw; } + } + + Triangles(Int4 nbvx,Geometry & G) :Gh(G),BTh(*this){ + try { GeomToTriangles0(nbvx);} + catch(...) { this->~Triangles(); throw; } + } + + Triangles(Triangles &,Geometry * pGh=0,Triangles* pBTh=0,Int4 nbvxx=0 ); // COPY OPERATEUR + Triangles(const Triangles &,const int *flag,const int *bb); // truncature void SetIntCoor(const char * from =0); @@ -746,5 +746,5 @@ void ForceBoundary(); void Heap(); - void FindSubDomain(int ); + void FindSubDomain(int OutSide=0); Int4 ConsRefTriangle(Int4 *) const; void ShowHistogram() const; @@ -820,5 +820,4 @@ class Geometry { public: - int OnDisk; Int4 NbRef; // counter of ref on the this class if 0 we can delete @@ -832,5 +831,5 @@ Int4 NbOfCurves; int empty(){return (nbv ==0) && (nbt==0) && (nbe==0) && (NbSubDomains==0); } - GeometricalVertex * vertices; // data of vertices des sommets + GeometricalVertex* vertices; Triangle * triangles; GeometricalEdge * edges; @@ -859,5 +858,5 @@ Geometry() {EmptyGeometry();}// empty Geometry void AfterRead(); - Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();OnDisk=1;ReadGeometry(bamggeom,bamgopts);AfterRead();} + Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts) {EmptyGeometry();ReadGeometry(bamggeom,bamgopts);AfterRead();} const GeometricalVertex & operator[] (Int4 i) const { return vertices[i];}; @@ -1143,5 +1142,5 @@ inline void Adj(GeometricalEdge * & on,int &i) - {int j=i;i=on->SensAdj[i];on=on->Adj[j];} + {int j=i;i=on->DirAdj[i];on=on->Adj[j];} inline Real4 qualite(const Vertex &va,const Vertex &vb,const Vertex &vc) Index: /issm/trunk/src/c/Bamgx/QuadTree.h =================================================================== --- /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2876) +++ /issm/trunk/src/c/Bamgx/QuadTree.h (revision 2877) @@ -9,5 +9,4 @@ const IntQuad MaxISize = ( 1L << MaxDeep); - class Triangles; class Vertex; @@ -15,8 +14,6 @@ class QuadTree { public: - class QuadTreeBox { public: - long n; // if n < 4 => Vertex else => QuadTreeBox; union { @@ -24,8 +21,5 @@ Vertex * v[4]; }; - - - }; // end class QuadTreeBox ///////////////// - + }; // end class class StorageQuadTreeBox { public: @@ -41,28 +35,20 @@ return len*sizeof(QuadTreeBox)+sizeof(StorageQuadTreeBox)+ (n?n->SizeOf():0); } - }; // end class StorageQuadTreeBox - + }; // end class StorageQuadTreeBox * sb; - - long lenStorageQuadTreeBox; public: - QuadTreeBox * root; - Triangles *th; + QuadTreeBox* root; + Triangles* th; long NbQuadTreeBox,NbVertices; long NbQuadTreeBoxSearch,NbVerticesSearch; - Vertex * NearestVertex(Icoor1 i,Icoor1 j); - Vertex * NearestVertexWithNormal(Icoor1 i,Icoor1 j); // new version - Vertex * ToClose(Vertex & ,Real8 ,Icoor1,Icoor1); + Vertex* NearestVertex(Icoor1 i,Icoor1 j); + Vertex* NearestVertexWithNormal(Icoor1 i,Icoor1 j); + Vertex* ToClose(Vertex & ,Real8 ,Icoor1,Icoor1); long SizeOf() const {return sizeof(QuadTree)+sb->SizeOf();} - - void Add( Vertex & w); - QuadTreeBox* NewQuadTreeBox(){ - if(! (sb->bcbe)) - sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); - + if(! (sb->bcbe)) sb=new StorageQuadTreeBox(lenStorageQuadTreeBox,sb); if (!sb || (sb->bc->n != 0)){ throw ErrorException(__FUNCT__,exprintf("!sb || (sb->bc->n != 0)")); Index: /issm/trunk/src/c/Bamgx/objects/Geometry.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2876) +++ /issm/trunk/src/c/Bamgx/objects/Geometry.cpp (revision 2877) @@ -314,5 +314,5 @@ edges[i].tg[0]=zero2; edges[i].tg[1]=zero2; - edges[i].SensAdj[0] = edges[i].SensAdj[1] = -1; + edges[i].DirAdj[0] = edges[i].DirAdj[1] = -1; edges[i].Adj[0] = edges[i].Adj[1] = 0; edges[i].flag = 0; @@ -483,206 +483,31 @@ /*Methods*/ - /*FUNCTION Geometry::EmptyGeometry(){{{1*/ - void Geometry::EmptyGeometry() { - OnDisk=0; - NbRef=0; - name =0; - quadtree=0; - curves=0; - // edgescomponante=0; - triangles=0; - edges=0; - vertices=0; - NbSubDomains=0; - // nbtf=0; - // BeginOfCurve=0; - nbiv=nbv=nbvx=0; - nbe=nbt=nbtx=0; - NbOfCurves=0; - // BeginOfCurve=0; - subdomains=0; - MaximalAngleOfCorner = 10*Pi/180; - } - /*}}}1*/ - /*FUNCTION Geometry::Contening{{{1*/ - GeometricalEdge* Geometry::Contening(const R2 P, GeometricalEdge * start) const { - GeometricalEdge* on =start,* pon=0; - // walk with the cos on geometry - int k=0; - while(pon != on){ - k++; - pon = on; - if (k>=100){ - throw ErrorException(__FUNCT__,exprintf("k>=100")); - } - R2 A= (*on)[0]; - R2 B= (*on)[1]; - R2 AB = B-A; - R2 AP = P-A; - R2 BP = P-B; - if ( (AB,AP) < 0) - on = on->Adj[0]; - else if ( (AB,BP) > 0) - on = on->Adj[1]; - else - return on; - } - return on; - } - /*}}}1*/ - /*FUNCTION Geometry::Geometry::ProjectOnCurve {{{1*/ - GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const { - Real8 save_s=s; - int NbTry=0; -retry: - s=save_s; - GeometricalEdge* on = e.on; - if (!on){ - throw ErrorException(__FUNCT__,exprintf("!on")); - } - if (!e[0].on || !e[1].on){ - throw ErrorException(__FUNCT__,exprintf("!e[0].on || !e[1].on")); - } - const Vertex &v0=e[0],&v1=e[1]; - V.m = Metric(1.0-s, v0,s, v1); - const int mxe =100; - GeometricalEdge *ge[mxe+1]; - int sensge[mxe+1]; - Real8 lge[mxe+1]; - int bge=mxe/2,tge=bge; - ge[bge] = e.on; - sensge[bge]=1; - - R2 V0 = v0,V1=v1,V01=V1-V0; - VertexOnGeom vg0= *v0.on, vg1=*v1.on; - - // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; - GeometricalEdge * eg0=on, *eg1=on; - R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; - int OppositeSens = (V01,AB) < 0; - int sens0=0,sens1=1; - if (OppositeSens) - s=1-s,Exchange(vg0,vg1),Exchange(V0,V1); - while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0){ - if (bge<=0) { - // int kkk; - if(NbTry) { - printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n"); - printf("That bug might come from:\n"); - printf(" 1) a mesh edge contening more than %i geometrical edges\n",mxe/2); - printf(" 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before\n"); - printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n"); - throw ErrorException(__FUNCT__,exprintf("see above")); - } - NbTry++; - goto retry;} - GeometricalEdge* tmpge = eg0; - ge[--bge] =eg0 = eg0->Adj[sens0]; - if (bge<0 || bge>mxe){ - throw ErrorException(__FUNCT__,exprintf("bge<0 || bge>mxe")); - } - sens0 = 1-( sensge[bge] = tmpge->SensAdj[sens0]); - } - while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(sens1) != (GeometricalVertex*) vg1) { - if(tge>=mxe ) { - printf("WARNING: on the class triangles before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n"); - NbTry++; - if (NbTry<2) goto retry; - printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n"); - printf("That bug might come from:\n"); - printf(" 1) a mesh edge contening more than %i geometrical edges\n",mxe/2); - printf(" 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before\n"); - printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n"); - throw ErrorException(__FUNCT__,exprintf("see above")); - } - - GeometricalEdge* tmpge = eg1; - ge[++tge] =eg1 = eg1->Adj[sens1]; - sensge[tge]= sens1 = 1-tmpge->SensAdj[sens1]; - if (tge<0 || tge>mxe){ - throw ErrorException(__FUNCT__,exprintf("(tge<0 || tge>mxe)")); - } - } - - - if ( (*eg0)(sens0) == (GeometricalVertex*) vg0 ) - vg0 = VertexOnGeom( *(Vertex *) vg0,*eg0,sens0); - - if ( (*eg1)(sens1) == (GeometricalVertex*) vg1) - vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1); - - Real8 sg; - if (eg0 == eg1) { - register Real8 s0= vg0,s1=vg1; - sg = s0 * (1.0-s) + s * s1; - on=eg0;} - else { - R2 AA=V0,BB; - Real8 s0,s1; - int i=bge; - Real8 ll=0; - for(i=bge;imxe){ - throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); - } - BB = (*ge[i])[sensge[i]]; - lge[i]=ll += Norme2(AA-BB); - AA=BB ;} - lge[tge]=ll+=Norme2(AA-V1); - // search the geometrical edge - if (s>1.0){ - throw ErrorException(__FUNCT__,exprintf("s>1.0")); - } - Real8 ls= s*ll; - on =0; - s0 = vg0; - s1= sensge[bge]; - Real8 l0=0,l1; - i=bge; - while ( (l1=lge[i]) < ls ) { - if (i<0 || i>mxe){ - throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); - } - i++,s0=1-(s1=sensge[i]),l0=l1;} - on=ge[i]; - if (i==tge) - s1=vg1; - - s=(ls-l0)/(l1-l0); - sg = s0 * (1.0-s) + s * s1; - } - if (!on){ - throw ErrorException(__FUNCT__,exprintf("!on")); - } - V.r= on->F(sg); - GV=VertexOnGeom(V,*on,sg); - return on; - } - /*}}}1*/ /*FUNCTION Geometry::AfterRead(){{{1*/ - void Geometry::AfterRead() { + void Geometry::AfterRead(){ long int verbosity=0; - Int4 i,k=0; ; - int jj; // jj in [0,1] - Int4 * hv = new Int4 [ nbv]; - Int4 * ev = new Int4 [ 2 * nbe ]; - float * eangle = new float[ nbe ]; - double eps = 1e-20; + Int4 i,j,k; + int jj; + Int4* hv=new Int4[nbv]; + Int4* ev=new Int4[2*nbe]; + float* eangle=new float[nbe]; + double eps=1e-20; QuadTree quadtree; // to find same vertices - Vertex * v0 = vertices; - GeometricalVertex * v0g = (GeometricalVertex *) (void *) v0; - - for (i=0;ir - vertices[i]) < eps ) - { // link v & vertices[i] - // vieille ruse pour recuperer j - GeometricalVertex * vg = (GeometricalVertex *) (void *) v; - int j = vg-v0g; + Vertex* v0=vertices; + GeometricalVertex* v0g=(GeometricalVertex*) (void*)v0; + + k=0; + //link all vertices to themselves by default + for (i=0;i0)) + for (i=0;ir - vertices[i]) < eps ){ + // mama's old trick to get j + GeometricalVertex* vg = (GeometricalVertex*) (void*)v; + j=vg-v0g; if ( v != &(Vertex &) vertices[j]){ throw ErrorException(__FUNCT__,exprintf(" v != &(Vertex &) vertices[j]")); @@ -690,7 +515,7 @@ vertices[i].link = vertices + j; k++; - } + } else quadtree.Add(vertices[i]); - } + } if (k) { printf("number of distinct vertices= %i, over %i\n",nbv - k,nbv); @@ -702,5 +527,5 @@ } - // verification of cracked edge + // verification of cracked edge (to be completed) for (i=0;ir - edges[i].v[0]->r; + R2 v10=edges[i].v[1]->r - edges[i].v[0]->r; Real8 lv10 = Norme2(v10); if(lv10 == 0) { @@ -731,9 +557,8 @@ } eangle[i] = atan2(v10.y,v10.x) ; // angle in [ -Pi,Pi ] - for (jj=0;jj<2;jj++) - { // generation of list - Int4 v = Number(edges[i].v[jj]); - ev[k] = hv[v]; - hv[v] = k++; + for (jj=0;jj<2;jj++){ + Int4 v=Number(edges[i].v[jj]); + ev[k]=hv[v]; + hv[v]=k++; } } @@ -804,5 +629,5 @@ } edges[i1].Adj[j1] = edges + i; - edges[i1].SensAdj[j1] = jj; + edges[i1].DirAdj[j1] = jj; } @@ -819,5 +644,5 @@ if( ! edges[i].v[jj]->Corner()) { // not a Corner tg = edges[i].v[1-jj]->r - - edges[i].Adj[jj]->v[1-edges[i].SensAdj[jj]]->r; + - edges[i].Adj[jj]->v[1-edges[i].DirAdj[jj]]->r; ltg = Norme2(tg); tg = tg *(lAB/ltg),ltg=lAB; @@ -838,8 +663,6 @@ } // for (i=0;iRequired() ) break; - k0 = e->SensAdj[k1];// vertex in next edge + k0 = e->DirAdj[k1];// vertex in next edge e = e->Adj[k1]; // next edge @@ -887,10 +710,9 @@ } - if(step==0) { - curves = new Curve[NbOfCurves]; - } - } - for(int i=0;ilink; //GeometricalEdge * ee=curves[i].ee, *eqee=be->link; @@ -909,6 +731,6 @@ if(be->ReverseEqui()) curves[i].Reverse(); - } - } + } + } delete []ev; @@ -918,4 +740,180 @@ } /*}}}1*/ + /*FUNCTION Geometry::Contening{{{1*/ + GeometricalEdge* Geometry::Contening(const R2 P, GeometricalEdge * start) const { + GeometricalEdge* on =start,* pon=0; + // walk with the cos on geometry + int k=0; + while(pon != on){ + k++; + pon = on; + if (k>=100){ + throw ErrorException(__FUNCT__,exprintf("k>=100")); + } + R2 A= (*on)[0]; + R2 B= (*on)[1]; + R2 AB = B-A; + R2 AP = P-A; + R2 BP = P-B; + if ( (AB,AP) < 0) + on = on->Adj[0]; + else if ( (AB,BP) > 0) + on = on->Adj[1]; + else + return on; + } + return on; + } + /*}}}1*/ + /*FUNCTION Geometry::EmptyGeometry(){{{1*/ + void Geometry::EmptyGeometry() { + NbRef=0; + name =0; + quadtree=0; + curves=0; + // edgescomponante=0; + triangles=0; + edges=0; + vertices=0; + NbSubDomains=0; + // nbtf=0; + // BeginOfCurve=0; + nbiv=nbv=nbvx=0; + nbe=nbt=nbtx=0; + NbOfCurves=0; + // BeginOfCurve=0; + subdomains=0; + MaximalAngleOfCorner = 10*Pi/180; + } + /*}}}1*/ + /*FUNCTION Geometry::Geometry::ProjectOnCurve {{{1*/ + GeometricalEdge* Geometry::ProjectOnCurve(const Edge & e,Real8 s,Vertex &V,VertexOnGeom &GV ) const { + Real8 save_s=s; + int NbTry=0; +retry: + s=save_s; + GeometricalEdge* on = e.on; + if (!on){ + throw ErrorException(__FUNCT__,exprintf("!on")); + } + if (!e[0].on || !e[1].on){ + throw ErrorException(__FUNCT__,exprintf("!e[0].on || !e[1].on")); + } + const Vertex &v0=e[0],&v1=e[1]; + V.m = Metric(1.0-s, v0,s, v1); + const int mxe =100; + GeometricalEdge *ge[mxe+1]; + int sensge[mxe+1]; + Real8 lge[mxe+1]; + int bge=mxe/2,tge=bge; + ge[bge] = e.on; + sensge[bge]=1; + + R2 V0 = v0,V1=v1,V01=V1-V0; + VertexOnGeom vg0= *v0.on, vg1=*v1.on; + + // GeometricalEdge * eg0 = e.on,* eg1 = e.on, *eg=NULL; + GeometricalEdge * eg0=on, *eg1=on; + R2 Ag=(R2) (*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag; + int OppositeSens = (V01,AB) < 0; + int sens0=0,sens1=1; + if (OppositeSens) + s=1-s,Exchange(vg0,vg1),Exchange(V0,V1); + while (eg0 != (GeometricalEdge*) vg0 && (*eg0)(sens0) != (GeometricalVertex*) vg0){ + if (bge<=0) { + // int kkk; + if(NbTry) { + printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n"); + printf("That bug might come from:\n"); + printf(" 1) a mesh edge contening more than %i geometrical edges\n",mxe/2); + printf(" 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before\n"); + printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n"); + throw ErrorException(__FUNCT__,exprintf("see above")); + } + NbTry++; + goto retry;} + GeometricalEdge* tmpge = eg0; + ge[--bge] =eg0 = eg0->Adj[sens0]; + if (bge<0 || bge>mxe){ + throw ErrorException(__FUNCT__,exprintf("bge<0 || bge>mxe")); + } + sens0 = 1-( sensge[bge] = tmpge->DirAdj[sens0]); + } + while (eg1 != (GeometricalEdge*) vg1 && (*eg1)(sens1) != (GeometricalVertex*) vg1) { + if(tge>=mxe ) { + printf("WARNING: on the class triangles before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)\n"); + NbTry++; + if (NbTry<2) goto retry; + printf("Fatal Error: on the class triangles before call Geometry::ProjectOnCurve\n"); + printf("That bug might come from:\n"); + printf(" 1) a mesh edge contening more than %i geometrical edges\n",mxe/2); + printf(" 2) code bug : be sure that we call Triangles::SetVertexFieldOn() before\n"); + printf("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)\n"); + throw ErrorException(__FUNCT__,exprintf("see above")); + } + + GeometricalEdge* tmpge = eg1; + ge[++tge] =eg1 = eg1->Adj[sens1]; + sensge[tge]= sens1 = 1-tmpge->DirAdj[sens1]; + if (tge<0 || tge>mxe){ + throw ErrorException(__FUNCT__,exprintf("(tge<0 || tge>mxe)")); + } + } + + + if ( (*eg0)(sens0) == (GeometricalVertex*) vg0 ) + vg0 = VertexOnGeom( *(Vertex *) vg0,*eg0,sens0); + + if ( (*eg1)(sens1) == (GeometricalVertex*) vg1) + vg1 = VertexOnGeom( *(Vertex *) vg1,*eg1,sens1); + + Real8 sg; + if (eg0 == eg1) { + register Real8 s0= vg0,s1=vg1; + sg = s0 * (1.0-s) + s * s1; + on=eg0;} + else { + R2 AA=V0,BB; + Real8 s0,s1; + int i=bge; + Real8 ll=0; + for(i=bge;imxe){ + throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); + } + BB = (*ge[i])[sensge[i]]; + lge[i]=ll += Norme2(AA-BB); + AA=BB ;} + lge[tge]=ll+=Norme2(AA-V1); + // search the geometrical edge + if (s>1.0){ + throw ErrorException(__FUNCT__,exprintf("s>1.0")); + } + Real8 ls= s*ll; + on =0; + s0 = vg0; + s1= sensge[bge]; + Real8 l0=0,l1; + i=bge; + while ( (l1=lge[i]) < ls ) { + if (i<0 || i>mxe){ + throw ErrorException(__FUNCT__,exprintf("i<0 || i>mxe")); + } + i++,s0=1-(s1=sensge[i]),l0=l1;} + on=ge[i]; + if (i==tge) + s1=vg1; + + s=(ls-l0)/(l1-l0); + sg = s0 * (1.0-s) + s * s1; + } + if (!on){ + throw ErrorException(__FUNCT__,exprintf("!on")); + } + V.r= on->F(sg); + GV=VertexOnGeom(V,*on,sg); + return on; + } + /*}}}1*/ } Index: /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2876) +++ /issm/trunk/src/c/Bamgx/objects/QuadTree.cpp (revision 2877) @@ -51,9 +51,8 @@ NbVertices(0), NbQuadTreeBoxSearch(0), - NbVerticesSearch(0) - { - sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); - root=NewQuadTreeBox(); - } + NbVerticesSearch(0){ + sb =new StorageQuadTreeBox(lenStorageQuadTreeBox); + root=NewQuadTreeBox(); + } /*}}}1*/ /*FUNCTION QuadTree::~QuadTree(){{{1*/ @@ -65,9 +64,52 @@ /*Methods*/ + /*FUNCTION QuadTree::Add{{{1*/ + void QuadTree::Add( Vertex & w){ + QuadTreeBox ** pb , *b; + register long i=w.i.x, j=w.i.y,l=MaxISize; + pb = &root; + while( (b=*pb) && (b->n<0)){ + b->n--; + l >>= 1; + pb = &b->b[IJ(i,j,l)]; + } + if (b) { + if (b->n > 3 && b->v[3] == &w) return; + if (b->n > 2 && b->v[2] == &w) return; + if (b->n > 1 && b->v[1] == &w) return; + if (b->n > 0 && b->v[0] == &w) return; + } + if (l==0){ + throw ErrorException(__FUNCT__,exprintf("l==0")); + } + while ((b= *pb) && (b->n == 4)){ // the QuadTreeBox is full + Vertex *v4[4]; // copy of the QuadTreeBox vertices + + v4[0]= b->v[0]; + v4[1]= b->v[1]; + v4[2]= b->v[2]; + v4[3]= b->v[3]; + b->n = -b->n; // mark is pointer QuadTreeBox + b->b[0]=b->b[1]=b->b[2]=b->b[3]=0; // set empty QuadTreeBox ptr + l >>= 1; // div the size by 2 + for (register int k=0;k<4;k++){ // for the 4 vertices find the sub QuadTreeBox ij + register int ij; + register QuadTreeBox * bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,l)]; + if (!bb) + bb=b->b[ij]=NewQuadTreeBox(); // alloc the QuadTreeBox + bb->v[bb->n++] = v4[k]; + } + pb = &b->b[IJ(i,j,l)]; + } + if (!(b = *pb)) b=*pb= NewQuadTreeBox(); // alloc the QuadTreeBox + b->v[b->n++]=&w; // we add the vertex + NbVertices++; + } + /*}}}1*/ /*FUNCTION QuadTree::NearestVertex{{{1*/ Vertex* QuadTree::NearestVertex(Icoor1 i,Icoor1 j) { - QuadTreeBox * pb[ MaxDeep ]; - int pi[ MaxDeep ]; - Icoor1 ii[ MaxDeep ], jj [ MaxDeep]; + QuadTreeBox* pb[MaxDeep]; + int pi[MaxDeep]; + Icoor1 ii[MaxDeep], jj [MaxDeep]; register int l=0; // level register QuadTreeBox * b; @@ -80,12 +122,11 @@ Vertex *vn=0; - // init for optimisation --- + // init for optimization b = root; - register Int4 n0; - if (!root->n) - return vn; // empty tree - - while( (n0 = b->n) < 0) - { + register Int4 n0; + + if (!root->n) return vn; // empty tree + + while( (n0 = b->n) < 0){ // search the non empty // QuadTreeBox containing the point (i,j) @@ -93,6 +134,5 @@ register int k = IJ(iplus,jplus,hb2);// QuadTreeBox number of size hb2 contening i;j register QuadTreeBox * b0= b->b[k]; - if ( ( b0 == 0) || (b0->n == 0) ) - break; // null box or empty => break + if ( ( b0 == 0) || (b0->n == 0) ) break; // null box or empty => break NbQuadTreeBoxSearch++; b=b0; @@ -100,11 +140,8 @@ j0 += J_IJ(k,hb2); // j orign of QuadTreeBox hb = hb2; - } - - - if ( n0 > 0) - { - for(register int k=0;k 0){ + for(register int k=0;kv[k]->i; h0 = NORM(iplus,i2.x,jplus,i2.y); @@ -113,8 +150,9 @@ vn = b->v[k];} NbVerticesSearch++; - } + } return vn; - } - // general case ----- + } + + // general case pb[0]= b; pi[0]=b->n>0 ?(int) b->n : 4 ; @@ -124,10 +162,8 @@ do { b= pb[l]; - while (pi[l]--) - { + while (pi[l]--){ register int k = pi[l]; - if (b->n>0) // Vertex QuadTreeBox none empty - { + if (b->n>0){ // Vertex QuadTreeBox none empty NbVerticesSearch++; I2 i2 = b->v[k]->i; @@ -138,17 +174,14 @@ vn = b->v[k]; } - } - else // Pointer QuadTreeBox - { + } + else{ // Pointer QuadTreeBox register QuadTreeBox *b0=b; NbQuadTreeBoxSearch++; - if ((b=b->b[k])) - { + if ((b=b->b[k])){ hb >>=1 ; // div by 2 register Icoor1 iii = ii[l]+I_IJ(k,hb); register Icoor1 jjj = jj[l]+J_IJ(k,hb); - if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)) - { + if (INTER_SEG(iii,iii+hb,iplus-h,iplus+h) && INTER_SEG(jjj,jjj+hb,jplus-h,jplus+h)){ pb[++l]= b; pi[l]= b->n>0 ?(int) b->n : 4 ; @@ -156,15 +189,14 @@ jj[l]= jjj; - } - else - b=b0, hb <<=1 ; - } - else - b=b0; - } - } + } + else{ + b=b0, hb <<=1 ; + } + } + else b=b0; + } + } hb <<= 1; // mul by 2 } while (l--); - return vn; } @@ -280,4 +312,18 @@ } /*}}}1*/ + /*FUNCTION QuadTree::StorageQuadTreeBox::StorageQuadTreeBox{{{1*/ + QuadTree::StorageQuadTreeBox::StorageQuadTreeBox(long ll,StorageQuadTreeBox *nn) { + len = ll; + n = nn; + b = new QuadTreeBox[ll]; + for (int i = 0; i n<0)) - { - b->n--; - l >>= 1; - pb = &b->b[IJ(i,j,l)]; - } - if (b) { - if (b->n > 3 && b->v[3] == &w) return; - if (b->n > 2 && b->v[2] == &w) return; - if (b->n > 1 && b->v[1] == &w) return; - if (b->n > 0 && b->v[0] == &w) return; - } - if (l==0){ - throw ErrorException(__FUNCT__,exprintf("l==0")); - } - while ((b= *pb) && (b->n == 4)) // the QuadTreeBox is full - { - Vertex *v4[4]; // copy of the QuadTreeBox vertices - - v4[0]= b->v[0]; - v4[1]= b->v[1]; - v4[2]= b->v[2]; - v4[3]= b->v[3]; - b->n = -b->n; // mark is pointer QuadTreeBox - b->b[0]=b->b[1]=b->b[2]=b->b[3]=0; // set empty QuadTreeBox ptr - l >>= 1; // div the size by 2 - for (register int k=0;k<4;k++) // for the 4 vertices find the sub QuadTreeBox ij - { - register int ij; - register QuadTreeBox * bb = b->b[ij=IJ(v4[k]->i.x,v4[k]->i.y,l)]; - if (!bb) - bb=b->b[ij]=NewQuadTreeBox(); // alloc the QuadTreeBox - bb->v[bb->n++] = v4[k]; - } - pb = &b->b[IJ(i,j,l)]; - } - if (!(b = *pb)) - b=*pb= NewQuadTreeBox(); // alloc the QuadTreeBox - b->v[b->n++]=&w; // we add the vertex - NbVertices++; - } - /*}}}1*/ - /*FUNCTION QuadTree::StorageQuadTreeBox::StorageQuadTreeBox{{{1*/ - QuadTree::StorageQuadTreeBox::StorageQuadTreeBox(long ll,StorageQuadTreeBox *nn) { - len = ll; - n = nn; - b = new QuadTreeBox[ll]; - for (int i = 0; i det; + + if (( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) ){ + throw ErrorException(__FUNCT__,exprintf("infv=%g det=%g")); + } + + // if det3 do not exist then constuct det3 + if (!det3) { + det3 = det3local; // alloc + if ( infv<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 + det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ; + det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ; + det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}} + + + if (!det3[0]) izerodet=0,nbd0++; + if (!det3[1]) izerodet=1,nbd0++; + if (!det3[2]) izerodet=2,nbd0++; + + if (nbd0 >0 ) // point s on a egde or on a vertex + if (nbd0 == 1) { + iedge = OppositeEdge[izerodet]; + TriangleAdjacent 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 ) { + // add in outside triangle + Add(s,( Triangle *) ta); + return;} + }} + else { + printf("bug (%i): Bug double points in\n",nbd0); + throw ErrorException(__FUNCT__,exprintf("See above")); + } + + // remove de MarkUnSwap edge + t->SetUnMarkUnSwap(0); + t->SetUnMarkUnSwap(1); + t->SetUnMarkUnSwap(2); + + tt[0]= t; + tt[1]= &triangles[nbt++]; + tt[2]= &triangles[nbt++]; + + if (nbt>nbtx) { + throw ErrorException(__FUNCT__,exprintf("Not ebough triangles")); + } + + *tt[1]= *tt[2]= *t; + // gestion of the link + tt[0]->link=tt[1]; + tt[1]->link=tt[2]; + + (* tt[0])(OppositeVertex[0])=&s; + (* tt[1])(OppositeVertex[1])=&s; + (* tt[2])(OppositeVertex[2])=&s; + + tt[0]->det=det3[0]; + tt[1]->det=det3[1]; + tt[2]->det=det3[2]; + + // update adj des triangles externe + tt[0]->SetAdjAdj(0); + tt[1]->SetAdjAdj(1); + tt[2]->SetAdjAdj(2); + // update des adj des 3 triangle interne + const int i0 = 0; + const int i1= NextEdge[i0]; + const int i2 = PreviousEdge[i0]; + + tt[i0]->SetAdj2(i2,tt[i2],i0); + tt[i1]->SetAdj2(i0,tt[i0],i1); + tt[i2]->SetAdj2(i1,tt[i1],i2); + + tt[0]->SetTriangleContainingTheVertex(); + tt[1]->SetTriangleContainingTheVertex(); + tt[2]->SetTriangleContainingTheVertex(); + + + // swap if the point s is on a edge + if(izerodet>=0) { + int rswap =tt[izerodet]->swap(iedge); + + if (!rswap) { + throw ErrorException(__FUNCT__,exprintf("swap the point s is on a edge")); + } + } + } + /*}}}1*/ + /*FUNCTION Triangles::BoundAnisotropy{{{1*/ + void Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso) { + long int verbosity=0; + + double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100)); + if (verbosity > 1) printf(" BoundAnisotropy by %g\n",anisomax); + Real8 h1=1.e30,h2=1e-30,rx=0; + Real8 coef = 1./(anisomax*anisomax); + Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; + for (Int4 i=0;i2){ + printf(" input: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5),pow(h1,-0.5),pow(rx,0.5)); + printf(" output: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); + } + } + /*}}}1*/ /*FUNCTION Triangles::ConsGeometry{{{1*/ void Triangles::ConsGeometry(Real8 cutoffradian,int *equiedges) // construct a geometry if no geo @@ -772,5 +927,5 @@ if (nbe != edge4->nb()){ throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i, nbe4=%i",nbe,edge4->nb())); - } + } for (i=0;inb()+nbv); } - // check the consistant of edge[].adj and the geometrical required vertex + // check the consistant of edge[].adj and the geometrical required vertex + k=0; + kk=0; + Int4 it; + + for (i=0;i=0) // edge alone + kk++; + + k += kk; + kk=0; + if (k) { + // construction of the edges + nbe = k; + Edge * edgessave = edges; + edges = new Edge[nbe]; + k =0; + // construction of the edges + if(verbosity>4) printf(" Construction of the edges %i\n",nbe); + + for (i=0;i=0) // edge alone + { + it = st[i]/3; + j = (int) (st[i]%3); + add=k++; + } + + if (add>=0 && add < nbe) + { + + edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]]; + edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]]; + edges[add].on=0; + if (iref(),edges[add].v[1]->ref()); // no a good choice + } + } + if (k!=nbe){ + throw ErrorException(__FUNCT__,exprintf("k!=nbe")); + } + if (edgessave) delete [] edgessave; + } + + // construction of edges[].adj + for (i=0;icolor++; + + for (i=0;icolor,j0; + if(i0<0) + edges[i ].adj[ j ]=0; // Add FH Jan 2008 + if(i0==-1) + v->color=i*2+j; + else if (i0>=0) {// i and i0 edge are adjacent by the vertex v + j0 = i0%2; + i0 = i0/2; + if (v!=edges[i0 ].v[j0]){ + throw ErrorException(__FUNCT__,exprintf("v!=edges[i0 ].v[j0]")); + } + edges[i ].adj[ j ] =edges +i0; + edges[i0].adj[ j0] =edges +i ; + v->color = -3;} + } + // now reconstruct the sub domain info + if (NbSubDomains){ + throw ErrorException(__FUNCT__,exprintf("NbSubDomains should be 0")); + } + NbSubDomains=0; + + { + Int4 it; + // find all the sub domain + Int4 *colorT = new Int4[nbt]; + Triangle *tt,*t; + Int4 k; + for ( it=0;it0) + if( ( j=st[level]++) <3) + { + t = &triangles[st[level-1]]; + tt=t->TriangleAdj((int)j); + + if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) + { + colorT[jt]=NbSubDomains; + st[++level]=jt; + st[++level]=0; + k++; + } + } + else + level-=2; + NbSubDomains++; + } + } + if (verbosity> 3) printf(" The Number of sub domain = %i\n",NbSubDomains); + + Int4 isd; + subdomains = new SubDomain[NbSubDomains]; + for (isd=0;isd=0) // edge alone - kk++; - - k += kk; - kk=0; - if (k) { - // construction of the edges - nbe = k; - Edge * edgessave = edges; - edges = new Edge[nbe]; - k =0; - // construction of the edges - if(verbosity>4) printf(" Construction of the edges %i\n",nbe); - - for (i=0;i=0) // edge alone - { - it = st[i]/3; - j = (int) (st[i]%3); - add=k++; - } - - if (add>=0 && add < nbe) - { - - edges[add].v[0] = &triangles[it][VerticesOfTriangularEdge[j][0]]; - edges[add].v[1] = &triangles[it][VerticesOfTriangularEdge[j][1]]; - edges[add].on=0; - if (iref(),edges[add].v[1]->ref()); // no a good choice - } - } - if (k!=nbe){ - throw ErrorException(__FUNCT__,exprintf("k!=nbe")); + for (it=0;itcolor != triangles[it].color) && !subdomains[isd=colorT[it]].head) + { + subdomains[isd].head = triangles+it; + subdomains[isd].ref = triangles[it].color; + subdomains[isd].sens = j; // hack + subdomains[isd].edge = 0; + k++; + } + } + if (k!= NbSubDomains){ + throw ErrorException(__FUNCT__,exprintf("k!= NbSubDomains")); + } + + delete [] colorT; + + + } + delete [] st; + // now make the geometry + // 1 compress the vertices + Int4 * colorV = new Int4[nbv]; + for (i=0;i3) printf(" number of vertices = %i\n number of edges = %i\n",Gh.nbv,Gh.nbe); + NbVerticesOnGeomVertex = Gh.nbv; + VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; + NbVerticesOnGeomEdge =0; + VerticesOnGeomEdge =0; + { + Int4 j; + for (i=0;i=0) + { + + Vertex & v = Gh.vertices[j]; + v = vertices[i]; + v.color =0; + VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]); } - if (edgessave) delete [] edgessave; + + } + edge4= new SetOfEdges4(nbe,nbv); + + Real4 * len = new Real4[Gh.nbv]; + for(i=0;i=0) { + if( sens==0) + Gh.edges[i].SetEqui(); + else + Gh.edges[i].SetReverseEqui(); + Gh.edges[i].link= & Gh.edges[j]; + } + + } + if(requis) { // correction fevr 2009 JYU ... + Gh.edges[i].v[0]->SetRequired(); + Gh.edges[i].v[1]->SetRequired(); + Gh.edges[i].SetRequired(); // fin modif ... + } + R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r; + Real8 l12=Norme2(x12); + hmin = Min(hmin,l12); + + Gh.vertices[j1].color++; + Gh.vertices[j0].color++; + + len[j0]+= l12; + len[j1] += l12; + hmin = Min(hmin,l12); + Gh.edges[i].ref = edges[i].ref; + + k = edge4->addtrie(i0,i1); + if (k != i){ + throw ErrorException(__FUNCT__,exprintf("k != i")); + } + + } + + + for (i=0;i 0) + Gh.vertices[i].m= Metric(len[i] /(Real4) Gh.vertices[i].color); + else + Gh.vertices[i].m= Metric(hmin); + delete [] len; + for (i=0;ifindtrie(i0,i1); + if(k>=0) + { + subdomains[i].sens = (vertices + i0 == edges[k].v[0]) ? 1 : -1; + subdomains[i].edge = edges+k; + Gh.subdomains[i].edge = Gh.edges + k; + Gh.subdomains[i].sens = subdomains[i].sens; + Gh.subdomains[i].ref = subdomains[i].ref; } - - // construction of edges[].adj - for (i=0;icolor++; - - for (i=0;icolor,j0; - if(i0<0) - edges[i ].adj[ j ]=0; // Add FH Jan 2008 - if(i0==-1) - v->color=i*2+j; - else if (i0>=0) {// i and i0 edge are adjacent by the vertex v - j0 = i0%2; - i0 = i0/2; - if (v!=edges[i0 ].v[j0]){ - throw ErrorException(__FUNCT__,exprintf("v!=edges[i0 ].v[j0]")); - } - edges[i ].adj[ j ] =edges +i0; - edges[i0].adj[ j0] =edges +i ; - v->color = -3;} - } - // now reconstruct the sub domain info - if (NbSubDomains){ - throw ErrorException(__FUNCT__,exprintf("NbSubDomains should be 0")); - } - NbSubDomains=0; - - { - Int4 it; - // find all the sub domain - Int4 *colorT = new Int4[nbt]; - Triangle *tt,*t; - Int4 k; - for ( it=0;it0) - if( ( j=st[level]++) <3) - { - t = &triangles[st[level-1]]; - tt=t->TriangleAdj((int)j); - - if ( ! t->Locked(j) && tt && (colorT[jt = Number(tt)] == -1) && ( tt->color==kolor)) - { - colorT[jt]=NbSubDomains; - st[++level]=jt; - st[++level]=0; - k++; - } - } - else - level-=2; - NbSubDomains++; - } - } - if (verbosity> 3) printf(" The Number of sub domain = %i\n",NbSubDomains); - - Int4 isd; - subdomains = new SubDomain[NbSubDomains]; - for (isd=0;isdcolor != triangles[it].color) && !subdomains[isd=colorT[it]].head) - { - subdomains[isd].head = triangles+it; - subdomains[isd].ref = triangles[it].color; - subdomains[isd].sens = j; // hack - subdomains[isd].edge = 0; - k++; - } - } - if (k!= NbSubDomains){ - throw ErrorException(__FUNCT__,exprintf("k!= NbSubDomains")); - } - - delete [] colorT; - - - } - delete [] st; - // now make the geometry - // 1 compress the vertices - Int4 * colorV = new Int4[nbv]; - for (i=0;i3) printf(" number of vertices = %i\n number of edges = %i\n",Gh.nbv,Gh.nbe); - NbVerticesOnGeomVertex = Gh.nbv; - VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; - NbVerticesOnGeomEdge =0; - VerticesOnGeomEdge =0; - { - Int4 j; - for (i=0;i=0) - { - - Vertex & v = Gh.vertices[j]; - v = vertices[i]; - v.color =0; - VerticesOnGeomVertex[j] = VertexOnGeom(vertices[i], Gh.vertices[j]); - } - - } - edge4= new SetOfEdges4(nbe,nbv); - - Real4 * len = new Real4[Gh.nbv]; - for(i=0;i=0) { - if( sens==0) - Gh.edges[i].SetEqui(); - else - Gh.edges[i].SetReverseEqui(); - Gh.edges[i].link= & Gh.edges[j]; - } - - } - if(requis) { // correction fevr 2009 JYU ... - Gh.edges[i].v[0]->SetRequired(); - Gh.edges[i].v[1]->SetRequired(); - Gh.edges[i].SetRequired(); // fin modif ... - } - R2 x12 = Gh.vertices[j0].r-Gh.vertices[j1].r; - Real8 l12=Norme2(x12); - hmin = Min(hmin,l12); - - Gh.vertices[j1].color++; - Gh.vertices[j0].color++; - - len[j0]+= l12; - len[j1] += l12; - hmin = Min(hmin,l12); - Gh.edges[i].ref = edges[i].ref; - - k = edge4->addtrie(i0,i1); - if (k != i){ - throw ErrorException(__FUNCT__,exprintf("k != i")); - } - - } - - - for (i=0;i 0) - Gh.vertices[i].m= Metric(len[i] /(Real4) Gh.vertices[i].color); - else - Gh.vertices[i].m= Metric(hmin); - delete [] len; - for (i=0;ifindtrie(i0,i1); - if(k>=0) - { - subdomains[i].sens = (vertices + i0 == edges[k].v[0]) ? 1 : -1; - subdomains[i].edge = edges+k; - Gh.subdomains[i].edge = Gh.edges + k; - Gh.subdomains[i].sens = subdomains[i].sens; - Gh.subdomains[i].ref = subdomains[i].ref; - } - else - throw ErrorException(__FUNCT__,exprintf("%i should be >=0")); - } - - delete edge4; - delete [] colorV; - // -- unset adj - for (i=0;i=0")); + } + + delete edge4; + delete [] colorV; + // -- unset adj + for (i=0;ibe; - tA=vA.onbe->abcisse; - } - else { - throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve On Vertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA))); - } - - if (vB.vint == IsVertexOnVertex){ - pB=vB.onbv; - } - else if(vB.vint == IsVertexOnEdge){ - pB=vB.onbe->be; - tB=vB.onbe->abcisse; - } - else { - throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve On Vertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB))); - } - Edge * e = &BhAB; - if (!pA || !pB || !e){ - throw ErrorException(__FUNCT__,exprintf("!pA || !pB || !e")); - } - // be carefull the back ground edge e is on same geom edge - // of the initiale edge def by the 2 vertex A B; - //check Is a background Mesh; - if (e=BTh.edges+BTh.nbe){ - throw ErrorException(__FUNCT__,exprintf("e=BTh.edges+BTh.nbe")); - } - // walk on BTh edge - //not finish ProjectOnCurve with BackGround Mesh); - // 1 first find a back ground edge contening the vertex A - // 2 walk n back gound boundary to find the final vertex B - - if( vA.vint == IsVertexOnEdge) - { // find the start edge - e = vA.onbe->be; - - } - else if (vB.vint == IsVertexOnEdge) - { - theta = 1-theta; - Exchange(tA,tB); - Exchange(pA,pB); - Exchange(pvA,pvB); - Exchange(A,B); - e = vB.onbe->be; - - } - else{ // do the search by walking - throw ErrorException(__FUNCT__,exprintf("case not supported yet")); - } - - // find the direction of walking with sens of edge and pA,PB; - R2 AB=B-A; - - Real8 cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB); - int kkk=0; - int sens = (cosE01AB>0) ? 1 : 0; - - // Real8 l=0; // length of the edge AB - Real8 abscisse = -1; - - for (int cas=0;cas<2;cas++){ - // 2 times algo: - // 1 for computing the length l - // 2 for find the vertex - int iii; - Vertex *v0=pvA,*v1; - Edge *neee,*eee; - Real8 lg =0; // length of the curve - Real8 te0; - // we suppose take the curve's abcisse - for ( eee=e,iii=sens,te0=tA; - eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ; - neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) { - - kkk=kkk+1; - if (kkk>=100){ - throw ErrorException(__FUNCT__,exprintf("kkk>=100")); + /*FUNCTION Triangles::ConsRefTriangle{{{1*/ + Int4 Triangles::ConsRefTriangle(Int4* reft) const { + long int verbosity=0; + register Triangle *t0,*t; + register Int4 k=0, num; + for (Int4 it=0;it=nbt){ + throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt")); } - if (!eee){ - throw ErrorException(__FUNCT__,exprintf("!eee")); + reft[num]=i; + } + while (t0 != (t=t->link)); + } + return k; + } + /*}}}1*/ + /*FUNCTION Triangles::Crack{{{1*/ + int Triangles::Crack() { + if (NbCrackedEdges!=0 && NbCrackedVertices<=0);{ + throw ErrorException(__FUNCT__,exprintf("NbCrackedEdges!=0 && NbCrackedVertices<=0")); + } + for (int i=0;iCracked()) k++; + if( k==0) return 0; + CurrentTh = this; + printf(" number of Cracked Edges = %i\n",k); + NbCrackedEdges =k; + CrackedEdges = new CrackedEdge[k]; + // new edge + Edge * e = new Edge[ nbe + k]; + + // copy + for (i=0;iCracked()) + { + e[nbe] = e[i]; + // return the edge + e[nbe].v[0] = e[i].v[1]; + e[nbe].v[1] = e[i].v[0]; + e[nbe].on = e[i].on->link ; // fqux + CrackedEdges[k++]=CrackedEdge(edges,i,nbe); + nbe++; } - Real8 lg0 = lg; - Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); - lg += dp; - if (cas && abscisse <= lg) { // ok we find the geom edge - Real8 sss = (abscisse-lg0)/dp; - Real8 thetab = te0*(1-sss)+ sss*iii; - if (thetab<0 || thetab>1){ - throw ErrorException(__FUNCT__,exprintf("thetab<0 || thetab>1")); + ReMakeTriangleContainingTheVertex() ; + // + int nbcrakev =0; + Vertex *vlast = vertices + nbv; + Vertex *vend = vertices + nbvx; // end of array + for (int iv=0;iv=3)){ + throw ErrorException(__FUNCT__,exprintf("!tbegin || (i<0) || (i>=3)")); + } + // turn around the vertex v + TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge + int k=0; + do { + int kv = VerticesOfTriangularEdge[ta][1]; + k++; + Triangle * tt (ta); + if ( ta.Cracked() ) + { + TriangleAdjacent tta=(ta.Adj()); + if (!tta.Cracked()){ + throw ErrorException(__FUNCT__,exprintf("!tta.Cracked()")); + } + if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart + if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use + kk++;// number of cracked edge view + } + if ( tt->link ) { // if good triangles store the value + int it = Number(tt); + if (it>=nt){ + throw ErrorException(__FUNCT__,exprintf("(it>=nt)")); + } + (*tt)(kv)= vv; // Change the vertex of triangle + if(vvReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref + // tt->SetTriangleContainingTheVertex(); + kkk++; + } else if (kk) { // crack + boundary + if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use } - BR = VertexOnEdge(&R,eee,thetab); - return Gh.ProjectOnCurve(*eee,thetab,R,GR); - } + + ta = Next(ta).Adj(); + } while ( (tbegin != ta)); + if (!k){ + throw ErrorException(__FUNCT__,exprintf("!k")); + } + if (kc) nbcrakev++; } - // we find the end - if (v1 != pvB){ - if (( void*) v1 == pB) - tB = iii; - - Real8 lg0 = lg; - if (!eee){ - throw ErrorException(__FUNCT__,exprintf("!eee")); - } - v1 = pvB; - Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); - lg += dp; - abscisse = lg*theta; - if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end - { // ok we find the geom edge - Real8 sss = (abscisse-lg0)/dp; - Real8 thetab = te0*(1-sss)+ sss*tB; - if (thetab<0 || thetab>1){ - throw ErrorException(__FUNCT__,exprintf("thetab<0 || thetab>1")); - } - BR = VertexOnEdge(&R,eee,thetab); - return Gh.ProjectOnCurve(*eee,thetab,R,GR); - } - } - abscisse = lg*theta; - - } - throw ErrorException(__FUNCT__,exprintf("Big bug...")); - return 0; // just for the compiler - } - /*}}}1*/ - /*FUNCTION Triangles::MakeQuadrangles{{{1*/ - void Triangles::MakeQuadrangles(double costheta){ - long int verbosity=0; - - if (verbosity>2) printf("MakeQuadrangles costheta = %g\n",costheta); - - if (costheta >1) { - if (verbosity>5) printf(" do nothing: costheta > 1\n"); - } - - - Int4 nbqq = (nbt*3)/2; - DoubleAndInt4 *qq = new DoubleAndInt4[nbqq]; - - Int4 i,ij; - int j; - Int4 k=0; - for (i=0;i=costheta) - qq[k++].i3j=i*3+j; - // sort qq - HeapSort(qq,k); - - Int4 kk=0; - for (ij=0;ij=costheta) - triangles[i].SetHidden(j),kk++; - } - NbOfQuad = kk; - if (verbosity>2){ - printf(" number of quadrilaterals = %i\n",NbOfQuad); - printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); - printf(" number of outside triangles = %i\n",NbOutT); - } - delete [] qq; - } - /*}}}1*/ - /*FUNCTION Triangles::SplitElement{{{1*/ - int Triangles::SplitElement(int choice){ - long int verbosity=0; - - Direction NoDirOfSearch; - const int withBackground = &BTh != this && &BTh; - - ReNumberingTheTriangleBySubDomain(); - //int nswap =0; - const Int4 nfortria( choice ? 4 : 6); - if(withBackground) - { - BTh.SetVertexFieldOn(); - SetVertexFieldOnBTh(); - } - else - BTh.SetVertexFieldOn(); - - Int4 newnbt=0,newnbv=0; - Int4 * kedge = 0; - Int4 newNbOfQuad=NbOfQuad; - Int4 * ksplit= 0, * ksplitarray=0; - Int4 kkk=0; - int ret =0; - if (nbvxaddtrie(Number(edges[i][0]),Number(edges[i][1]))); - if (ong) // a geometrical edges - { - if (withBackground){ - // walk on back ground mesh - // newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); - // a faire -- difficile - // the first PB is to now a background edge between the 2 vertices - if (!edgesGtoB){ - throw ErrorException(__FUNCT__,exprintf("!edgesGtoB")); - } - ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].on)], - edges[i][0],edges[i][1],0.5,vertices[k], - newVertexOnBThEdge[kvb], - newVerticesOnGeomEdge[kvg++]); - vertices[k].ReferenceNumber= edges[i].ref; - vertices[k].DirOfSearch = NoDirOfSearch; - ; - // get the Info on background mesh - Real8 s = newVertexOnBThEdge[kvb]; - Vertex & bv0 = newVertexOnBThEdge[kvb][0]; - Vertex & bv1 = newVertexOnBThEdge[kvb][1]; - // compute the metrix of the new points - vertices[k].m = Metric(1-s,bv0,s,bv1); - kvb++; - } - else - { - ong=Gh.ProjectOnCurve(edges[i], - 0.5,vertices[k],newVerticesOnGeomEdge[kvg++]); - // vertices[k].i = toI2( vertices[k].r); - vertices[k].ReferenceNumber = edges[i].ref; - vertices[k].DirOfSearch = NoDirOfSearch; - vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); - } - } - else // straigth line edge --- - { - vertices[k].r = ((R2) edges[i][0] + (R2) edges[i][1] )*0.5; - vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); - vertices[k].on = 0; - } - //vertices[k].i = toI2( vertices[k].r); - R2 AB = vertices[k].r; - R2 AA = (A+AB)*0.5; - R2 BB = (AB+B)*0.5; - vertices[k].ReferenceNumber = edges[i].ref; - vertices[k].DirOfSearch = NoDirOfSearch; - - newedges[ie].on = Gh.Contening(AA,ong); - newedges[ie++].v[1]=vertices+k; - - newedges[ie]=edges[i]; - newedges[ie].adj[0]=newedges + ie -1; - newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ; - newedges[ie].on = Gh.Contening(BB,ong); - newedges[ie++].v[0]=vertices+k; - k++; - } - if (edgesGtoB) delete [] edgesGtoB; - edgesGtoB=0; - - newnbv=k; - newNbVerticesOnGeomEdge=kvg; - if (newnbv> nbvx) goto Error;// bug - - nbv = k; - - - kedge = new Int4[3*nbt+1]; - ksplitarray = new Int4[nbt+1]; - ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0] - - for (i=0;i<3*nbt;i++) - kedge[i]=-1; - - // - - for (i=0;i= lastT) - t.SetAdj2(j,0,0);// unset adj - const Vertex & v0 = t[VerticesOfTriangularEdge[j][0]]; - const Vertex & v1 = t[VerticesOfTriangularEdge[j][1]]; - Int4 ke =edge4->findtrie(Number(v0),Number(v1)); - if (ke>0) - { - Int4 ii = Number(tt); - int jj = ta; - Int4 ks = ke + nbvold; - kedge[3*i+j] = ks; - if (ii=0){ - // warning PB roundoff error - if (t.link && ( (aa=Area2( A.r , t[1].r , t[2].r )) < 0.0 - || (bb=Area2( t[0].r , A.r , t[2].r )) < 0.0 - || (cc=Area2( t[0].r , t[1].r , A.r )) < 0.0)){ - printf("%i not in triangle %i In= %i %i %i %i %i\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd); - throw ErrorException(__FUNCT__,exprintf("Number of triangles with P2 interpolation Problem")); + + if ( nbcrakev ) + for (int iec =0;iec < NbCrackedEdges; iec ++) + CrackedEdges[iec].Set(); + + // set the ref + NbCrackedVertices = nbcrakev; + // int nbvo = nbv; + nbv = vlast - vertices; + int nbnewv = nbv - nbv; // nb of new vrtices + if (nbcrakev && verbosity > 1 ) printf(" number of Cracked vertices = %i, number of created vertices = %i\n",nbcrakev,nbnewv); + // all the new vertices are on geometry + if (nbnewv) + { // + Int4 n = nbnewv+NbVerticesOnGeomVertex; + Int4 i,j,k; + VertexOnGeom * vog = new VertexOnGeom[n]; + for ( i =0; i= LastOld) + { // a new vertex + Int4 old = v->ReferenceNumber ; // the old same vertex + Int4 i = ( v - LastOld); + // if the old is on vertex => warning + // else the old is on edge => ok + vog[i] = vog[old]; + // vog[i].mv = v; + //g[i].ge = ; + //og[i].abcisse = ; + } + } } - else { - if (tt.link && ( (aa=Area2( A.r , tt[1].r , tt[2].r )) < 0 - || (bb=Area2( tt[0].r , A.r , tt[2].r )) < 0 - || (cc=Area2( tt[0].r , tt[1].r , A.r )) < 0)){ - printf("%i not in triangle %i In= %i %i %i %i %i\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd); - throw ErrorException(__FUNCT__,exprintf("Number of triangles with P2 interpolation Problem")); - } - } - } + + NbVerticesOnGeomVertex = n; } - } - - for (i=0;ifindtrie(Number(v0),Number(v1)); - if (ke<0) // new - { - if (&tt) // internal triangles all the boundary - { // new internal edges - Int4 ii = Number(tt); - int jj = ta; - - kedge[3*i+j]=k;// save the vertex number - kedge[3*ii+jj]=k; - if (k=0 - kedge[3*i+j]=nbvold+ke; - kkk[nbsplitedge++]=j;// previously splited - } - } - else - kkk[nbsplitedge++]=j;// previously splited - - } - if (nbinvisible>=2){ - throw ErrorException(__FUNCT__,exprintf("nbinvisible>=2")); - } - switch (nbsplitedge) { - case 0: ksplit[i]=10; newnbt++; break; // nosplit - case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2 - case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3 - case 3: - if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4; - else ksplit[i]=10*nfortria,newnbt+=nfortria; - break; - } - if (ksplit[i]<40){ - throw ErrorException(__FUNCT__,exprintf("ksplit[i]<40")); - } - } - // now do the element split - newNbOfQuad = 4*NbOfQuad; - nbv = k; - kkk = nbt; - ksplit[-1] = nbt; - // look on old true triangles - - for (i=0;i=7 || kk<=0){ - throw ErrorException(__FUNCT__,exprintf("kk>=7 || kk<=0")); - } - - // def the numbering k (edge) i vertex - int k0 = ke; - int k1 = NextEdge[k0]; - int k2 = PreviousEdge[k0]; - int i0 = OppositeVertex[k0]; - int i1 = OppositeVertex[k1]; - int i2 = OppositeVertex[k2]; - - Triangle &t0=triangles[i]; - Vertex * v0=t0(i0); - Vertex * v1=t0(i1); - Vertex * v2=t0(i2); - - if (nbmkadj>=10){ - throw ErrorException(__FUNCT__,exprintf("nbmkadj>=10")); - } - // -------------------------- - TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2)); - // save the flag Hidden - int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)}; - // un set all adj -- save Hidden flag -- - t0.SetAdj2(0,0,hid[0]); - t0.SetAdj2(1,0,hid[1]); - t0.SetAdj2(2,0,hid[2]); - // -- remake - switch (kk) { - case 1: break;// nothing - case 2: // - { - Triangle &t1=triangles[kkk++]; - t1=t0; - if (kedge[3*i+i0]<0){ - throw ErrorException(__FUNCT__,exprintf("kedge[3*i+i0]<0")); - } - Vertex * v3 = vertices + kedge[3*i+k0]; - - t0(i2) = v3; - t1(i1) = v3; - t0.SetAllFlag(k2,0); - t1.SetAllFlag(k1,0); - } - break; - case 3: // - { - Triangle &t1=triangles[kkk++]; - Triangle &t2=triangles[kkk++]; - t2=t1=t0; - if (kedge[3*i+k1]<0){ - throw ErrorException(__FUNCT__,exprintf("kedge[3*i+k1]<0")); - } - if (kedge[3*i+k2]<0){ - throw ErrorException(__FUNCT__,exprintf("kedge[3*i+k2]<0")); - } - - Vertex * v01 = vertices + kedge[3*i+k2]; - Vertex * v02 = vertices + kedge[3*i+k1]; - t0(i1) = v01; - t0(i2) = v02; - t1(i2) = v02; - t1(i0) = v01; - t2(i0) = v02; - t0.SetAllFlag(k0,0); - t1.SetAllFlag(k1,0); - t1.SetAllFlag(k0,0); - t2.SetAllFlag(k2,0); - } - break; - case 4: // - case 6: // split in 4 - { - Triangle &t1=triangles[kkk++]; - Triangle &t2=triangles[kkk++]; - Triangle &t3=triangles[kkk++]; - t3=t2=t1=t0; - if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){ - throw ErrorException(__FUNCT__,exprintf("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0")); - } - Vertex * v12 = vertices + kedge[3*i+k0]; - Vertex * v02 = vertices + kedge[3*i+k1]; - Vertex * v01 = vertices + kedge[3*i+k2]; - t0(i1) = v01; - t0(i2) = v02; - t0.SetAllFlag(k0,hid[k0]); - - t1(i0) = v01; - t1(i2) = v12; - t0.SetAllFlag(k1,hid[k1]); - - t2(i0) = v02; - t2(i1) = v12; - t2.SetAllFlag(k2,hid[k2]); - - t3(i0) = v12; - t3(i1) = v02; - t3(i2) = v01; - - t3.SetAllFlag(0,hid[0]); - t3.SetAllFlag(1,hid[1]); - t3.SetAllFlag(2,hid[2]); - - if ( kk == 6) - { - - Triangle &t4=triangles[kkk++]; - Triangle &t5=triangles[kkk++]; - - t4 = t3; - t5 = t3; - - t0.SetHidden(k0); - t1.SetHidden(k1); - t2.SetHidden(k2); - t3.SetHidden(0); - t4.SetHidden(1); - t5.SetHidden(2); - - if (nbv < nbvx ) - { - vertices[nbv].r = ((R2) *v01 + (R2) *v12 + (R2) *v02 ) / 3.0; - vertices[nbv].ReferenceNumber =0; - vertices[nbv].DirOfSearch =NoDirOfSearch; - //vertices[nbv].i = toI2(vertices[nbv].r); - Real8 a3[]={1./3.,1./3.,1./3.}; - vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m); - Vertex * vc = vertices +nbv++; - t3(i0) = vc; - t4(i1) = vc; - t5(i2) = vc; - - } - else - goto Error; - } - - } - break; - } - - // save all the new triangles - mkadj[nbmkadj++]=i; - Int4 jj; - if (t0.link) - for (jj=nbt;jj13){// 13 = 6 + 4 + - throw ErrorException(__FUNCT__,exprintf("nbmkadj>13")); - } - - if (kk==6) newNbOfQuad+=3; - for (jj=ksplit[i-1];jj= nbtx) goto Error; // bug - if (nbv>= nbvx) goto Error; // bug - // generation of the new triangles - - SetIntCoor("In SplitElement"); - - ReMakeTriangleContainingTheVertex(); - if(withBackground) - BTh.ReMakeTriangleContainingTheVertex(); - - delete [] edges; - edges = newedges; - nbe = newnbe; - NbOfQuad = newNbOfQuad; - - for (i=0;i2){ - printf(" number of quadrilaterals = %i\n",NbOfQuad); - printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); - printf(" number of outside triangles = %i\n",NbOutT); - } - - CurrentTh=OCurrentTh; - return 0; //ok - -Error: - nbv = nbvold; - nbt = nbtold; - NbOutT = NbOutTold; - // cleaning memory --- - delete newedges; - if (ksplitarray) delete [] ksplitarray; - if (kedge) delete [] kedge; - if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge; - if (edge4) delete edge4; - if(newVertexOnBThEdge) delete [] newVertexOnBThEdge; - - - CurrentTh= OCurrentTh; - return ret; // ok - } - /*}}}1*/ - /*FUNCTION Triangles::swap{{{1*/ - int Triangle::swap(Int2 a,int koption){ - if(a/4 !=0) return 0;// arete lock or MarkUnSwap - - register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent - register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles - if(a2/4 !=0) return 0; // arete lock or MarkUnSwap - - register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]]; - register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]]; - register Vertex *s1=t1->ns[OppositeVertex[a1]]; - register Vertex *s2=t2->ns[OppositeVertex[a2]]; - - Icoor2 det1=t1->det , det2=t2->det ; - Icoor2 detT = det1+det2; - Icoor2 detA = Abs(det1) + Abs(det2); - Icoor2 detMin = Min(det1,det2); - - int OnSwap = 0; - // si 2 triangle infini (bord) => detT = -2; - if (sa == 0) {// les deux triangles sont frontieres - det2=bamg::det(s2->i,sb->i,s1->i); - OnSwap = det2 >0;} - else if (sb == 0) { // les deux triangles sont frontieres - det1=bamg::det(s1->i,sa->i,s2->i); - OnSwap = det1 >0;} - else if(( s1 != 0) && (s2 != 0) ) { - det1 = bamg::det(s1->i,sa->i,s2->i); - det2 = detT - det1; - OnSwap = (Abs(det1) + Abs(det2)) < detA; - - Icoor2 detMinNew=Min(det1,det2); - // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test - if (! OnSwap &&(detMinNew>0)) { - OnSwap = detMin ==0; - if (! OnSwap) { - int kopt = koption; - while (1) - if(kopt) { - // critere de Delaunay pure isotrope - register Icoor2 xb1 = sb->i.x - s1->i.x, - x21 = s2->i.x - s1->i.x, - yb1 = sb->i.y - s1->i.y, - y21 = s2->i.y - s1->i.y, - xba = sb->i.x - sa->i.x, - x2a = s2->i.x - sa->i.x, - yba = sb->i.y - sa->i.y, - y2a = s2->i.y - sa->i.y; - register double - cosb12 = double(xb1*x21 + yb1*y21), - cosba2 = double(xba*x2a + yba*y2a) , - sinb12 = double(det2), - sinba2 = double(t2->det); - - - // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) - OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); - break; - } - else - { - // critere de Delaunay anisotrope - Real8 som; - I2 AB=(I2) *sb - (I2) *sa; - I2 MAB2=((I2) *sb + (I2) *sa); - R2 MAB(MAB2.x*0.5,MAB2.y*0.5); - I2 A1=(I2) *s1 - (I2) *sa; - I2 D = (I2) * s1 - (I2) * sb ; - R2 S2(s2->i.x,s2->i.y); - R2 S1(s1->i.x,s1->i.y); - { - Metric M=s1->m; - R2 ABo = M.Orthogonal(AB); - R2 A1o = M.Orthogonal(A1); - // (A+B)+ x ABo = (S1+B)/2+ y A1 - // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 - double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); - double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 - if (Abs(d) > dd*1.e-3) { - R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); - som = M(C - S2)/M(C - S1); - } else - {kopt=1;continue;} - - } - { - Metric M=s2->m; - R2 ABo = M.Orthogonal(AB); - R2 A1o = M.Orthogonal(A1); - // (A+B)+ x ABo = (S1+B)/2+ y A1 - // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 - double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); - double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 - if(Abs(d) > dd*1.e-3) { - R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); - som += M(C - S2)/M(C - S1); - } else - {kopt=1;continue;} - } - OnSwap = som < 2; - break; - } - - } // OnSwap - } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) - } - if( OnSwap ) - bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); - else { - NbUnSwap ++; - t1->SetMarkUnSwap(a1); - } - return OnSwap; - } - /*}}}1*/ - /*FUNCTION Triangles::MetricAt{{{1*/ - Metric Triangles::MetricAt(const R2 & A) const { - I2 a = toI2(A); - Icoor2 deta[3]; - Triangle * t =FindTriangleContening(a,deta); - if (t->det <0) { // outside - double ba,bb; - TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ; - return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));} - else { // inside - Real8 aa[3]; - Real8 s = deta[0]+deta[1]+deta[2]; - aa[0]=deta[0]/s; - aa[1]=deta[1]/s; - aa[2]=deta[2]/s; - return Metric(aa,(*t)[0],(*t)[1],(*t)[2]); - } - } - /*}}}1*/ - /*FUNCTION Triangles::Add{{{1*/ - void Triangles::Add( Vertex & s,Triangle * t, Icoor2 * det3) { - // ------------------------------------------- - // s2 - // ! - // /|\ ! - // / | \ ! - // / | \ ! - // tt1 / | \ tt0 ! - // / |s \ ! - // / . \ ! - // / . ` \ ! - // / . ` \ ! - // ---------------- ! - // s0 tt2 s1 - //-------------------------------------------- - - Triangle * tt[3]; // the 3 new Triangles - Vertex &s0 = (* t)[0], &s1=(* t)[1], &s2=(* t)[2]; - Icoor2 det3local[3]; - int infv = &s0 ? (( &s1 ? ( &s2 ? -1 : 2) : 1 )) : 0; - // infv = ordre of the infini vertex (null) - register int nbd0 =0; // number of zero det3 - register int izerodet=-1,iedge; // izerodet = egde contening the vertex s - Icoor2 detOld = t->det; - - if (( infv <0 ) && (detOld <0) || ( infv >=0 ) && (detOld >0) ){ - throw ErrorException(__FUNCT__,exprintf("infv=%g det=%g")); - } - - // if det3 do not exist then constuct det3 - if (!det3) { - det3 = det3local; // alloc - if ( infv<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 - det3[0]= &s0 ? -1 : bamg::det(s ,s1,s2) ; - det3[1]= &s1 ? -1 : bamg::det(s0,s ,s2) ; - det3[2]= &s2 ? -1 : bamg::det(s0,s1,s ) ;}} - - - if (!det3[0]) izerodet=0,nbd0++; - if (!det3[1]) izerodet=1,nbd0++; - if (!det3[2]) izerodet=2,nbd0++; - - if (nbd0 >0 ) // point s on a egde or on a vertex - if (nbd0 == 1) { - iedge = OppositeEdge[izerodet]; - TriangleAdjacent 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 ) { - // add in outside triangle - Add(s,( Triangle *) ta); - return;} - }} - else { - printf("bug (%i): Bug double points in\n",nbd0); - throw ErrorException(__FUNCT__,exprintf("See above")); - } - - // remove de MarkUnSwap edge - t->SetUnMarkUnSwap(0); - t->SetUnMarkUnSwap(1); - t->SetUnMarkUnSwap(2); - - tt[0]= t; - tt[1]= &triangles[nbt++]; - tt[2]= &triangles[nbt++]; - - if (nbt>nbtx) { - throw ErrorException(__FUNCT__,exprintf("Not ebough triangles")); - } - - *tt[1]= *tt[2]= *t; - // gestion of the link - tt[0]->link=tt[1]; - tt[1]->link=tt[2]; - - (* tt[0])(OppositeVertex[0])=&s; - (* tt[1])(OppositeVertex[1])=&s; - (* tt[2])(OppositeVertex[2])=&s; - - tt[0]->det=det3[0]; - tt[1]->det=det3[1]; - tt[2]->det=det3[2]; - - // update adj des triangles externe - tt[0]->SetAdjAdj(0); - tt[1]->SetAdjAdj(1); - tt[2]->SetAdjAdj(2); - // update des adj des 3 triangle interne - const int i0 = 0; - const int i1= NextEdge[i0]; - const int i2 = PreviousEdge[i0]; - - tt[i0]->SetAdj2(i2,tt[i2],i0); - tt[i1]->SetAdj2(i0,tt[i0],i1); - tt[i2]->SetAdj2(i1,tt[i1],i2); - - tt[0]->SetTriangleContainingTheVertex(); - tt[1]->SetTriangleContainingTheVertex(); - tt[2]->SetTriangleContainingTheVertex(); - - - // swap if the point s is on a edge - if(izerodet>=0) { - int rswap =tt[izerodet]->swap(iedge); - - if (!rswap) { - throw ErrorException(__FUNCT__,exprintf("swap the point s is on a edge")); - } - } - } - /*}}}1*/ - /*FUNCTION Triangles::SplitInternalEdgeWithBorderVertices{{{1*/ - Int4 Triangles::SplitInternalEdgeWithBorderVertices(){ - Int4 NbSplitEdge=0; - SetVertexFieldOn(); - Int4 it; - Int4 nbvold=nbv; - long int verbosity=2; - for (it=0;itlink) { - printf("problem inserting point\n"); - } - - quadtree->Add(vi); - if (!tcvi || tcvi->det<0){// internal - throw ErrorException(__FUNCT__,exprintf("!tcvi || tcvi->det < 0")); - } - Add(vi,tcvi,dete); - NbSwap += vi.Optim(1); - iv++; - // } - } - if (verbosity>3) { - printf(" number of points: %i\n",iv); - printf(" number of swap to split internal edges with border vertices: %i\n",NbSwap); - nbv = iv; - } - } - if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold)); - if (verbosity>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge); - - return NbSplitEdge; - } - /*}}}1*/ - /*FUNCTION Triangles::InsertNewPoints{{{1*/ - Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) { - long int verbosity=2; - Real8 seuil= 1.414/2 ;// for two close point - Int4 i; - // insertion part --- - - const Int4 nbvnew = nbv-nbvold; - if (verbosity>5) printf(" Try to Insert the %i new points\n",nbvnew); - Int4 NbSwap=0; - Icoor2 dete[3]; - - // construction d'un ordre aleatoire - if (! nbvnew) - return 0; - if (nbvnew) { - const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ; - Int4 k3 = rand()%nbvnew ; - for (Int4 is3=0; is3ReferenceNumber=i; - } - // be carefull - Int4 iv = nbvold; - for (i=nbvold;iToClose(vi,seuil,hi,hj)) - { - // a good new point - Vertex & vj = vertices[iv]; - Int4 j = vj.ReferenceNumber; - if ( &vj!= ordre[j]){ - throw ErrorException(__FUNCT__,exprintf("&vj!= ordre[j]")); - } - if(i!=j) - { // for valgring - Exchange(vi,vj); - Exchange(ordre[j],ordre[i]); - } - vj.ReferenceNumber=0; - Triangle *tcvj = FindTriangleContening(vj.i,dete); - if (tcvj && !tcvj->link){ - throw ErrorException(__FUNCT__,exprintf("problem inserting point")); - } - quadtree->Add(vj); - Add(vj,tcvj,dete); - NbSwap += vj.Optim(1); - iv++; - } - } - if (verbosity>3) { - printf(" number of new points: %i\n",iv); - printf(" number of to close (?) points: %i\n",nbv-iv); - printf(" number of swap after: %i\n",NbSwap); - } - nbv = iv; - } - - for (i=nbvold;i3) printf(" NbSwap=%i\n",NbSwap); - - NbTSwap += NbSwap ; - return nbv-nbvold; - } - /*}}}1*/ - /*FUNCTION Triangles::NewPoints{{{1*/ - void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex) { - long int verbosity=2; - Int4 nbtold(nbt),nbvold(nbv); - if (verbosity>2) printf(" Triangles::NewPoints\n"); - Int4 i,k; - int j; - Int4 *first_np_or_next_t = new Int4[nbtx]; - Int4 NbTSwap =0; - // insert old point - nbtold = nbt; - - if (KeepBackVertex && (&Bh != this) && (nbv+Bh.nbv< nbvx)){ - // Bh.SetVertexFieldOn(); - for (i=0;i= nbt - // the list of test triangle is - // the next traingle on i is -first_np_or_next_t[i] - int iter=0; - // Big loop - do { - iter++; - nbtold = nbt; - nbvold = nbv; - - // default size of IntersectionTriangle - - i=Headt; - next_t=-first_np_or_next_t[i]; - for(t=&triangles[i];i=nbt ){ - throw ErrorException(__FUNCT__,exprintf("i<0 || i>=nbt")); - } - first_np_or_next_t[i] = iter; - for(j=0;j<3;j++) - { // for each edge - TriangleAdjacent tj(t,j); - Vertex & vA = * tj.EdgeVertex(0); - Vertex & vB = * tj.EdgeVertex(1); - - if (!t->link) continue;// boundary - if (t->det <0) continue; - if (t->Locked(j)) continue; - - TriangleAdjacent tadjj = t->Adj(j); - Triangle * ta= tadjj; - - if (ta->det <0) continue; - - R2 A = vA; - R2 B = vB; - - k=Number(ta); - - if(first_np_or_next_t[k]==iter) // this edge is done before - continue; // next edge of the triangle - - //const Int4 NbvOld = nbv; - lIntTria.SplitEdge(Bh,A,B); - lIntTria.NewPoints(vertices,nbv,nbvx); - } // end loop for each edge - - }// for triangle - - if (!InsertNewPoints(nbvold,NbTSwap)) - break; - - for (i=nbtold;it, EdgesVertexTriangle[s->vint][1]); - Triangle * tbegin= (Triangle*) ta; - Int4 kt; - do { - kt = Number((Triangle*) ta); - if (first_np_or_next_t[kt]>0) - first_np_or_next_t[kt]=-Headt,Headt=kt; - if (ta.EdgeVertex(0)!=s){ - throw ErrorException(__FUNCT__,exprintf("ta.EdgeVertex(0)!=s")); - } - ta = Next(Adj(ta)); - } while ( (tbegin != (Triangle*) ta)); - } - - } while (nbv!=nbvold); - - delete [] first_np_or_next_t; - - Int4 NbSwapf =0,NbSwp; - - // bofbof - - - NbSwp = NbSwapf; - for (i=0;i2) printf(" Insert initial %i vertices\n",nbv); - Triangles * OldCurrentTh =CurrentTh; - - CurrentTh=this; - SetIntCoor(); - Int4 i; - for (i=0;ii, ordre[1]->i, ordre[i]->i ) == 0;){ - if ( ++i >= nbv) { - throw ErrorException(__FUNCT__,exprintf("all the vertices are aligned")); - } - } - // echange i et 2 dans ordre afin - // que les 3 premiers ne soit pas aligne - Exchange( ordre[2], ordre[i]); - - // on ajoute un point a l'infini pour construire le maillage - // afin d'avoir une definition simple des aretes frontieres - nbt = 2; - - // on construit un maillage trivale forme - // d'une arete et de 2 triangles - // construit avec le 2 aretes orientes et - Vertex * v0=ordre[0], *v1=ordre[1]; - - triangles[0](0) = 0; // sommet pour infini - triangles[0](1) = v0; - triangles[0](2) = v1; - - triangles[1](0) = 0;// sommet pour infini - triangles[1](2) = v0; - triangles[1](1) = v1; - const int e0 = OppositeEdge[0]; - const int e1 = NextEdge[e0]; - const int e2 = PreviousEdge[e0]; - triangles[0].SetAdj2(e0, &triangles[1] ,e0); - triangles[0].SetAdj2(e1, &triangles[1] ,e2); - triangles[0].SetAdj2(e2, &triangles[1] ,e1); - - triangles[0].det = -1; // faux triangles - triangles[1].det = -1; // faux triangles - - triangles[0].SetTriangleContainingTheVertex(); - triangles[1].SetTriangleContainingTheVertex(); - - triangles[0].link=&triangles[1]; - triangles[1].link=&triangles[0]; - - // nbtf = 2; - if ( !quadtree ) quadtree = new QuadTree(this,0); - quadtree->Add(*v0); - quadtree->Add(*v1); - - //the vertices are added one by one - Int4 NbSwap=0; - - - if (verbosity>3) printf(" Begining of insertion process...\n"); - - for (Int4 icount=2; icounti,dete); - quadtree->Add(*vi); - Add(*vi,tcvi,dete); - NbSwap += vi->Optim(1,0); - } - if (verbosity>3) { - printf(" NbSwap of insertion: %i\n",NbSwap); - printf(" NbSwap/nbv: %i\n",NbSwap/nbv); - printf(" NbUnSwap: %i\n",NbUnSwap); - printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); - } - NbUnSwap = 0; - // construction d'un ordre aleatoire - // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ; -#ifdef NBLOOPOPTIM - - k3 = rand()%nbv ; - for (int is4=0; is4Optim(0,0); - if (verbosity>3) { - printf(" Optim Loop: %i\n",Nbloop); - printf(" NbSwap/nbv: %i\n",NbSwap/nbv); - printf(" NbUnSwap: %i\n",NbUnSwap); - printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); - } - NbUnSwap = 0; - if(!NbSwap) break; - } - ReMakeTriangleContainingTheVertex(); - // because we break the TriangleContainingTheVertex -#endif - CurrentTh=OldCurrentTh; - } - /*}}}1*/ + SetVertexFieldOn(); + + if (vlast >= vend) { + throw ErrorException(__FUNCT__,exprintf("Not enougth vertices: to crack the mesh we need %i vertices",nbv)); + } + CurrentTh = CurrentThOld; + return NbCrackedVertices; + } + /*}}}1*/ /*FUNCTION Triangles::ForceBoundary{{{1*/ void Triangles::ForceBoundary() { @@ -2497,6 +1495,297 @@ } /*}}}1*/ + /*FUNCTION Triangles::FillHoleInMesh{{{1*/ + void Triangles::FillHoleInMesh() { + + Triangles* OldCurrentTh =CurrentTh; + CurrentTh=this; + + int verbosity=0; + + // generation of the integer coordinate + { + + // find extrema coordinates of vertices pmin,pmax + Int4 i; + if(verbosity>2) printf(" Filling holes in mesh of %i vertices\n",nbv); + + //initialize ordre + if (!ordre){ + throw ErrorException(__FUNCT__,exprintf("!ordre")); + } + for (i=0;iaddtrie(Number(edges[i][0]),Number(edges[i][1]))); + } + if (kk != nbe) { + throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe)); + } + + // + for (i=0;iaddtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), + Number(triangles[i][VerticesOfTriangularEdge[j][1]])); + Int4 invisible = triangles[i].Hidden(j); + if(st[k]==-1){ + st[k]=3*i+j; + } + else if(st[k]>=0) { + if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){ + throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)))")); + } + + triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); + if (invisible) triangles[i].SetHidden(j); + if (k5) { + printf(" info on Mesh %s:\n",name); + printf(" - number of vertices = %i \n",nbv); + printf(" - number of triangles = %i \n",nbt); + printf(" - number of given edges = %i \n",nbe); + printf(" - number of all edges = %i \n" ,edge4->nb()); + printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); + } + + // check the consistant of edge[].adj and the geometrical required vertex + Int4 k=0; + for (i=0;inb();i++){ + if (st[i] >=0){ // edge alone + if (i < nbe) { + Int4 i0=edge4->i(i); + ordre[i0] = vertices+i0; + Int4 i1=edge4->j(i); + ordre[i1] = vertices+i1; + } + else { + k=k+1; + if (k <20) { + throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i))); + } + } + } + } + if(k != 0) { + throw ErrorException(__FUNCT__,exprintf("%i boundary edges are not defined as edges",k)); + } + + /* mesh generation with boundary points*/ + Int4 nbvb = 0; + for (i=0;ii, ordre[1]->i, ordre[i]->i ) == 0;) + if ( ++i >= nbvb) { + throw ErrorException(__FUNCT__,exprintf("FillHoleInMesh: All the vertices are aligned")); + } + Exchange( ordre[2], ordre[i]); + + Vertex * v0=ordre[0], *v1=ordre[1]; + + + triangles[0](0) = 0; // sommet pour infini + triangles[0](1) = v0; + triangles[0](2) = v1; + + triangles[1](0) = 0;// sommet pour infini + triangles[1](2) = v0; + triangles[1](1) = v1; + const int e0 = OppositeEdge[0]; + const int e1 = NextEdge[e0]; + const int e2 = PreviousEdge[e0]; + triangles[0].SetAdj2(e0, &triangles[1] ,e0); + triangles[0].SetAdj2(e1, &triangles[1] ,e2); + triangles[0].SetAdj2(e2, &triangles[1] ,e1); + + triangles[0].det = -1; // faux triangles + triangles[1].det = -1; // faux triangles + + triangles[0].SetTriangleContainingTheVertex(); + triangles[1].SetTriangleContainingTheVertex(); + + triangles[0].link=&triangles[1]; + triangles[1].link=&triangles[0]; + + if (!quadtree ) + delete quadtree; // ->ReInitialise(); + + quadtree = new QuadTree(this,0); + quadtree->Add(*v0); + quadtree->Add(*v1); + + // We add the vertices one by one + Int4 NbSwap=0; + for (Int4 icount=2; icounti,dete); + quadtree->Add(*vi); + Add(*vi,tcvi,dete); + NbSwap += vi->Optim(1,1); + } + + // inforce the boundary + TriangleAdjacent ta(0,0); + Int4 nbloss = 0,knbe=0; + for ( i = 0; i < nbe; i++){ + if (st[i] >=0){ // edge alone => on border ... FH oct 2009 + Vertex & a=edges[i][0], & b = edges[i][1]; + if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1. + { + knbe++; + if (ForceEdge(a,b,ta)<0) + nbloss++; + } + } + } + if(nbloss) { + throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe)); + } + + FindSubDomain(1); + // remove all the hole + // remove all the good sub domain + Int4 krm =0; + for (i=0;iaddtrie(v0?Number(v0):nbv,v1? Number(v1):nbv); + if (st[k]<0){ + throw ErrorException(__FUNCT__,exprintf("st[k]<0")); + } + tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); + ta.SetLock(); + st[k]=-2-st[k]; + } + } + } + } + Int4 NbTfillHoll =0; + for (i=0;isavenbtx ){ + throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx")); + } + // copy of the outside triangles in saveTriangles + for (i=0;i=0) { + savetriangles[savenbt]=triangles[i]; + savetriangles[savenbt].link=0; + savenbt++; + } + } + // gestion of the adj + k =0; + Triangle * tmax = triangles + nbt; + for (i=0;i= triangles && ta < tmax) + { + ta= savetriangles + ta->color; + ti.SetAdj2(j,ta,aa); + if(lck) ti.SetLocked(j); + } + } + } + // OutSidesTriangles = triangles; + // Int4 NbOutSidesTriangles = nbt; + + // restore triangles; + nbt=savenbt; + nbtx=savenbtx; + delete [] triangles; + delete [] subdomains; + triangles = savetriangles; + subdomains = savesubdomains; + if (k) { + throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k)); + } + FindSubDomain(); + + delete edge4; + delete [] st; + for (i=0;iAdd(vertices[i]); + + SetVertexFieldOn(); + + for (i=0;iIsRequiredVertex()) { + throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)")); + } + } + CurrentTh=OldCurrentTh; + } + /*}}}1*/ /*FUNCTION Triangles::FindSubDomain{{{1*/ - void Triangles::FindSubDomain(int OutSide=0) { + void Triangles::FindSubDomain(int OutSide) { long int verbosity=0; @@ -2747,162 +2036,347 @@ } /*}}}1*/ - /*FUNCTION Triangles::ReNumberingVertex{{{1*/ - void Triangles::ReNumberingVertex(Int4 * renu) { - // warning be carfull because pointeur - // from on mesh to over mesh - // -- so do ReNumbering a the beginning - Vertex * ve = vertices+nbv; - Int4 it,ie,i; - - for ( it=0;it=vertices && v < ve) - VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)]; - } - - for (i=0;i< NbVerticesOnGeomEdge;i++) - { - Vertex *v =VerticesOnGeomEdge[i].mv; - if (v>=vertices && v < ve) - VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)]; - } - - for (i=0;i< NbVertexOnBThVertex;i++) - { - Vertex *v=VertexOnBThVertex[i].v; - if (v>=vertices && v < ve) - VertexOnBThVertex[i].v=vertices+renu[Number(v)]; - } - - for (i=0;i< NbVertexOnBThEdge;i++) - { - Vertex *v=VertexOnBThEdge[i].v; - if (v>=vertices && v < ve) - VertexOnBThEdge[i].v=vertices+renu[Number(v)]; - } - - // move the Vertices without a copy of the array - // be carefull not trivial code - Int4 j; - for ( it=0;it= 0) // a new sub cycle - { - i=it; - Vertex ti=vertices[i],tj; - while ( (j=renu[i]) >= 0) - { // i is old, and j is new - renu[i] = -1-renu[i]; // mark - tj = vertices[j]; // save new - vertices[j]= ti; // new <- old - i=j; // next - ti = tj; - } - } - if (quadtree) - { delete quadtree; - quadtree = new QuadTree(this); - } - for ( it=0;itNearestVertex(B.x,B.y) ; + + if (! a || !a->t ) { + if (a) { + printf("TriangleConteningTheVertex vertex number %i, another call to ReMakeTriangleContainingTheVertex was required\n", Number(a)); + } + throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContening")); + } + if (a=vertices+nbv){ + throw ErrorException(__FUNCT__,exprintf("a=vertices+nbv")); + } + t = a->t; + if (t=triangles+nbt){ + throw ErrorException(__FUNCT__,exprintf("t=triangles+nbt")); + } + } + Icoor2 detop ; + int kkkk =0; // number of test triangle + + while ( t->det < 0){ // the initial triangles is outside + int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; + if (k0<0){ // k0 the NULL vertex + throw ErrorException(__FUNCT__,exprintf("k0<0")); + } + 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 + return t; + t = t->TriangleAdj(OppositeEdge[k0]); + kkkk++; + if (kkkk>=2){ + throw ErrorException(__FUNCT__,exprintf("kkkk>=2")); + } + } + + jj=0; + detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B); + + while(t->det > 0 ) { + kkkk++; + if (kkkk>=2000){ + throw ErrorException(__FUNCT__,exprintf("kkkk>=2000")); + } + j= OppositeVertex[jj]; + dete[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 + 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; + // 0 => ok + // 1 => go in way 1 + // 2 => two way go in way 1 or 2 randomly + + if (k==0) break; + if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]); + if ( k>=3){ + throw ErrorException(__FUNCT__,exprintf("k>=3")); + } + TriangleAdjacent t1 = t->Adj(jj=ii[0]); + if ((t1.det() < 0 ) && (k == 2)) + t1 = t->Adj(jj=ii[1]); + t=t1; + j=t1;// for optimisation we now the -det[OppositeVertex[j]]; + detop = -dete[OppositeVertex[jj]]; + jj = j; + } + + if (t->det<0) // outside triangle + dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop; + // NbOfTriangleSearchFind += kkkk; + return t; } /*}}}1*/ - /*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/ - void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) { - long int verbosity=0; - Int4 *renu= new Int4[nbt]; - register Triangle *t0,*t,*te=triangles+nbt; - register Int4 k=0,it,i,j; - - for ( it=0;it= nbt ){ - throw ErrorException(__FUNCT__,exprintf("kt<0 || kt >= nbt")); + /*FUNCTION Triangles::GeomToTriangles0{{{1*/ + void Triangles::GeomToTriangles0(Int4 inbvx) { + /*Generate mesh from geometry*/ + + Gh.NbRef++;// add a ref to GH + + Int4 i,NbOfCurves=0,NbNewPoints,NbEdgeCurve; + Real8 lcurve,lstep,s; + + R2 AB; + GeometricalVertex *a,*b; + Vertex *va,*vb; + GeometricalEdge * e; + + //initialize this + PreInit(inbvx); + + int background = &BTh != this; + nbv=0; + NbVerticesOnGeomVertex=0; + NbVerticesOnGeomEdge=0; + + //Compute bumber of vertices on geometrical vertex + for (i=0;i= nbvx) { + throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx)); + } + + //compute nbv (number of vertices) + nbv=0; + for (i=0;i Th + VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]); + nbv++; + } + } + + //check that edges is empty + if (edges){ + throw ErrorException(__FUNCT__,exprintf("edges")); + } + + // 2 step + // step=0 to compute the number of edges + alloc at end + // step=1 to construct the edges + for (int step=0;step<2;step++){ + + //initialize number of edges and number of edges max + Int4 nbex=0; + nbe=0; + Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; + + + Gh.UnMarkEdges(); + NbOfCurves = 0; + + //go through the edges of the geometry + for (i=0;iThe(); + b=ei(1)->The(); + if (!edges){ + throw ErrorException(__FUNCT__,exprintf("!edges")); + } + edges[nbe].v[0]=a->to; + edges[nbe].v[1]=b->to;; + edges[nbe].ref = e->ref; + edges[nbe].on = e; + edges[nbe].adj[0] = 0; + edges[nbe].adj[1] = 0; + nbe++; + } + } + } + else { // on curve ------ + for ( int kstep=0;kstep<= step;kstep++){ + // if 2nd step where 2 step + // 1 compute the length of the curve + // 2 create the points and edge + PreviousNewEdge=0; + NbNewPoints=0; + NbEdgeCurve=0; + if (nbvend>=nbvx){ + throw ErrorException(__FUNCT__,exprintf("nbvend>=nbvx")); + } + lcurve =0; + s = lstep; + int k=j; + e = & ei; + a=ei(k)->The(); + va = a->to; + e->SetMark(); + + // if SameGeo We have go in the background geometry + // to find the discretisation of the curve + for(;;){ + k = 1-k; + b= (*e)(k)->The(); + AB = b->r - a->r; + Metric MA = background ? BTh.MetricAt(a->r) :a->m ; + Metric MB = background ? BTh.MetricAt(b->r) :b->m ; + Real8 ledge = (MA(AB) + MB(AB))/2; + + const int MaxSubEdge = 10; + int NbSubEdge = 1; + Real8 lSubEdge[MaxSubEdge]; + R2 A,B; + if (ledge < 1.5) lSubEdge[0] = ledge; + else { + NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); + A= a->r; + Metric MAs =MA,MBs; + ledge = 0; + Real8 x =0, xstep= 1. / NbSubEdge; + for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { + x += xstep; + B = e->F(k ? x : 1-x); + MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB); + AB = A-B; + lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); + } + } + + Real8 lcurveb = lcurve+ ledge ; + while (lcurve<=s && s <= lcurveb && nbv < nbvend){ + // New points + + // Real8 aa=(lcurveb-s)/ledge; + // Real8 bb=(s-lcurve)/ledge; + + Real8 ss = s-lcurve; + // 1) find the SubEdge containing ss by dichotomie + int kk0=-1,kk1=NbSubEdge-1,kkk; + Real8 ll0=0,ll1=ledge,llk; + while (kk1-kk0>1) + { + if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) + kk1=kkk,ll1=llk; + else + kk0=kkk,ll0=llk;} + if (kk1 == kk0){ + throw ErrorException(__FUNCT__,exprintf("kk1 == kk0")); + } + + Real8 sbb = (ss-ll0 )/(ll1-ll0); + Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb; + + // new vertex on edge + vb = &vertices[nbv++]; + vb->m = Metric(aa,a->m,bb,b->m); + vb->ReferenceNumber = e->ref; + vb->DirOfSearch =NoDirOfSearch; + Real8 abcisse = k ? bb : aa; + vb->r = e->F( abcisse ); + VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); + + // to take into account the sens of the edge + s += lstep; + edges[nbe].v[0]=va; + edges[nbe].v[1]=vb; + edges[nbe].ref = e->ref; + edges[nbe].on = e; + edges[nbe].adj[0] = PreviousNewEdge; + if(PreviousNewEdge) + PreviousNewEdge->adj[1] = &edges[nbe]; + PreviousNewEdge = edges + nbe; + nbe++; + va = vb; + } + lcurve = lcurveb; + e->SetMark(); + a=b; + if (b->Required() ) break; + int kprev=k; + k = e->DirAdj[kprev];// next vertices + e = e->Adj[kprev]; + if (!e){ + throw ErrorException(__FUNCT__,exprintf("!e")); + } + }// for(;;) + vb = b->to; + NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1); + NbNewPoints = NbEdgeCurve-1; + if(!kstep){ + NbVerticesOnGeomEdge0 += NbNewPoints; + NbOfCurves++; + } + nbvend=nbv+NbNewPoints; + lstep = lcurve / NbEdgeCurve; + }// end of curve -- + if (edges) { // last edges of the curves + edges[nbe].v[0]=va; + edges[nbe].v[1]=vb; + edges[nbe].ref = e->ref; + edges[nbe].on = e; + edges[nbe].adj[0] = PreviousNewEdge; + edges[nbe].adj[1] = 0; + if(PreviousNewEdge) PreviousNewEdge->adj[1] = & edges[nbe]; + nbe++; + } + else nbe += NbEdgeCurve; + } // end on curve --- + } // if (edges[i][j].Corner()) + } } - if (renu[kt]!=-1){ - throw ErrorException(__FUNCT__,exprintf("renu[kt]!=-1")); + } // for (i=0;ilink)); - } - // take is same numbering if possible - if(justcompress) - for ( k=0,it=0;it=0 ) - renu[it]=k++; - - // put the outside triangles at the end - for ( it=0;it= 0) // a new sub cycle - { - i=it; - Triangle ti=triangles[i],tj; - while ( (j=renu[i]) >= 0) - { // i is old, and j is new - renu[i] = -1; // mark - tj = triangles[j]; // save new - triangles[j]= ti; // new <- old - i=j; // next - ti = tj; - } - } - delete [] renu; - nt = nbt - NbOutT; - - } - /*}}}1*/ - /*FUNCTION Triangles::ConsRefTriangle{{{1*/ - Int4 Triangles::ConsRefTriangle(Int4* reft) const { - long int verbosity=0; - register Triangle *t0,*t; - register Int4 k=0, num; - for (Int4 it=0;it=nbt){ - throw ErrorException(__FUNCT__,exprintf("num<0 || num>=nbt")); + if (VerticesOnGeomEdge){ + throw ErrorException(__FUNCT__,exprintf("VerticesOnGeomEdge")); } - reft[num]=i; - } - while (t0 != (t=t->link)); - } - return k; + edges = new Edge[nbex=nbe]; + if(NbVerticesOnGeomEdge0) + VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0]; + if (!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0){ + throw ErrorException(__FUNCT__,exprintf("!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0")); + } + // do the vertex on a geometrical vertex + NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge; + } + else if (NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0){ + throw ErrorException(__FUNCT__,exprintf("NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0")); + } + } + + Insert(); + ForceBoundary(); + FindSubDomain(); + + // NewPointsOld(*this) ; + NewPoints(*this,0) ; + CurrentTh = 0; } /*}}}1*/ @@ -3248,1802 +2722,2326 @@ } /*}}}1*/ - /*FUNCTION Triangles::GeomToTriangles0{{{1*/ - void Triangles::GeomToTriangles0(Int4 inbvx) { - Gh.NbRef++;// add a ref to GH - - - Int4 i,NbOfCurves=0,NbNewPoints,NbEdgeCurve; - Real8 lcurve, lstep,s; - - R2 AB; - GeometricalVertex *a,*b; - Vertex *va,*vb; - GeometricalEdge * e; - PreInit(inbvx); - int background = &BTh != this; - // int SameGeom = background && (&BTh.Gh == &Gh); - nbv = 0; - NbVerticesOnGeomVertex=0; - NbVerticesOnGeomEdge=0; - for (i=0;i= nbvx) { - throw ErrorException(__FUNCT__,exprintf("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,nbvx)); - } - for (i=0;i Th - VerticesOnGeomVertex[nbv]= VertexOnGeom(*Gh[i].to,Gh[i]); - nbv++; - } - - // Method in 2 step: 0 and 1 - // 1) compute de nb of edge - // 2) construct the edge - // generation of the curves - if (edges){ - throw ErrorException(__FUNCT__,exprintf("edges")); - } - // 2 step - // --step=0 to compute the number of edges + alloc at end - // --step=1 to construct the edges - for (int step=0;step<2;step++) - {// for (int step=0;step<2;step++) - Int4 nbex = 0; - nbe = 0; - Int4 NbVerticesOnGeomEdge0=NbVerticesOnGeomEdge; - Gh.UnMarkEdges(); - NbOfCurves = 0; - for (i=0;iThe(); - b=ei(1)->The(); - if (!edges){ - throw ErrorException(__FUNCT__,exprintf("!edges")); - } - edges[nbe].v[0]=a->to; - edges[nbe].v[1]=b->to;; - edges[nbe].ref = e->ref; - edges[nbe].on = e; - edges[nbe].adj[0] = 0; - edges[nbe].adj[1] = 0; - nbe++;} - } - else - { // on curve ------ - for ( int kstep=0;kstep<= step;kstep++) - { // begin for ( int kstep=0;kstep<= step;kstep++) - // if 2nd step where 2 step - // -- 1 compute le length of the curve - // -- create the points and edge - PreviousNewEdge=0; - NbNewPoints=0; - NbEdgeCurve=0; - if (nbvend>=nbvx){ - throw ErrorException(__FUNCT__,exprintf("nbvend>=nbvx")); - } - lcurve =0; - s = lstep; - int k=j; - e = & ei; - a=ei(k)->The(); - va = a->to; - e->SetMark(); - - // if SameGeo We have go in the background geometry - // to find the discretisation of the curve - - for(;;) - { - k = 1-k; - b= (*e)(k)->The(); - AB = b->r - a->r; - Metric MA = background ? BTh.MetricAt(a->r) :a->m ; - Metric MB = background ? BTh.MetricAt(b->r) :b->m ; - Real8 ledge = (MA(AB) + MB(AB))/2; - // - const int MaxSubEdge = 10; - int NbSubEdge = 1; - Real8 lSubEdge[MaxSubEdge]; - R2 A,B; - if (ledge < 1.5) - lSubEdge[0] = ledge; - else { - NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5)); - A= a->r; - Metric MAs =MA,MBs; - ledge = 0; - Real8 x =0, xstep= 1. / NbSubEdge; - for (int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) { - x += xstep; - B = e->F(k ? x : 1-x); - MBs= background ? BTh.MetricAt(B) :Metric(1-x, MA, x ,MB); - AB = A-B; - lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); - } - } - - Real8 lcurveb = lcurve+ ledge ; - while (lcurve<=s && s <= lcurveb && nbv < nbvend) - { - // New points - - // Real8 aa=(lcurveb-s)/ledge; - // Real8 bb=(s-lcurve)/ledge; - - Real8 ss = s-lcurve; - // 1) find the SubEdge containing ss by dichotomie - int kk0=-1,kk1=NbSubEdge-1,kkk; - Real8 ll0=0,ll1=ledge,llk; - while (kk1-kk0>1) - { - if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2])) - kk1=kkk,ll1=llk; - else - kk0=kkk,ll0=llk;} - if (kk1 == kk0){ - throw ErrorException(__FUNCT__,exprintf("kk1 == kk0")); - } - - Real8 sbb = (ss-ll0 )/(ll1-ll0); - Real8 bb = (kk1+sbb)/NbSubEdge, aa=1-bb; - - // new vertex on edge - vb = &vertices[nbv++]; - vb->m = Metric(aa,a->m,bb,b->m); - vb->ReferenceNumber = e->ref; - vb->DirOfSearch =NoDirOfSearch; - Real8 abcisse = k ? bb : aa; - vb->r = e->F( abcisse ); - VerticesOnGeomEdge[NbVerticesOnGeomEdge++]= VertexOnGeom(*vb,*e,abcisse); - - // to take in account the sens of the edge - - s += lstep; - edges[nbe].v[0]=va; - edges[nbe].v[1]=vb; - edges[nbe].ref = e->ref; - edges[nbe].on = e; - edges[nbe].adj[0] = PreviousNewEdge; - if(PreviousNewEdge) - PreviousNewEdge->adj[1] = &edges[nbe]; - PreviousNewEdge = edges + nbe; - nbe++; - va = vb; - } - lcurve = lcurveb; - e->SetMark(); - a=b; - if (b->Required() ) break; - int kprev=k; - k = e->SensAdj[kprev];// next vertices - e = e->Adj[kprev]; - if (!e){ - throw ErrorException(__FUNCT__,exprintf("!e")); - } - }// for(;;) - vb = b->to; - NbEdgeCurve = Max((Int4) (lcurve +0.5), (Int4) 1); - NbNewPoints = NbEdgeCurve-1; - if(!kstep) - { NbVerticesOnGeomEdge0 += NbNewPoints; - NbOfCurves++;} - - nbvend=nbv+NbNewPoints; - - lstep = lcurve / NbEdgeCurve; - } - // end of curve -- - if (edges) { // last edges of the curves - edges[nbe].v[0]=va; - edges[nbe].v[1]=vb; - edges[nbe].ref = e->ref; - edges[nbe].on = e; - edges[nbe].adj[0] = PreviousNewEdge; - edges[nbe].adj[1] = 0; - if(PreviousNewEdge) - PreviousNewEdge->adj[1] = & edges[nbe]; - - nbe++;} - else - nbe += NbEdgeCurve; - } // end on curve --- - } // if (edges[i][j].Corner()) - } // for (i=0;i H is computed with green formule + // otherwise => H is computed from P2 on 4T + const int dim = 2; + + long int verbosity=0; + + int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; + + // computation of the nb of field + Int4 ntmp = 0; + if (typsols) + { + for (Int4 i=0;i1e-30; + Real8 hmin = Max(hmin1,MinimalHmin()); + Real8 hmax = Min(hmax1,MaximalHmax()); + Real8 coef2 = 1/(coef*coef); + + if(verbosity>1) { + printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); + printf(" coef = %g\n",coef); + printf(" hmin = %g hmax = %g\n",hmin,hmax); + printf(" anisomax = %g nb Jacobi = %i, power = %i\n",anisomax,NbJacobi,power); + if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff); + else printf(" Absolute error\n"); } - /*}}}1*/ - /*FUNCTION Triangles::MakeGeometricalEdgeToEdge{{{1*/ - Edge** Triangles::MakeGeometricalEdgeToEdge() { - if (!Gh.nbe){ - throw ErrorException(__FUNCT__,exprintf("!Gh.nbe")); - } - Edge **e= new (Edge* [Gh.nbe]); - - Int4 i; - for ( i=0;ion; - e[Gh.Number(on)] = ei; + double *ss=(double*)s; + + double sA,sB,sC; + + Real8 *detT = new Real8[nbt]; + Real8 *Mmass= new Real8[nbv]; + Real8 *Mmassxx= new Real8[nbv]; + Real8 *dxdx= new Real8[nbv]; + Real8 *dxdy= new Real8[nbv]; + Real8 *dydy= new Real8[nbv]; + Real8 *workT= new Real8[nbt]; + Real8 *workV= new Real8[nbv]; + int *OnBoundary = new int[nbv]; + for (iv=0;ivlink) // no adj triangle => edge on boundary + OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1, + OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1, + nbb++; + } + + workT[i] = nbb; + Mmass[iA] += dett; + Mmass[iB] += dett; + Mmass[iC] += dett; + + if((nbb==0)|| !choice) + { + Mmassxx[iA] += dett; + Mmassxx[iB] += dett; + Mmassxx[iC] += dett; + } + } + else + workT[i]=-1; + + for (Int4 nusol=0;nusolon; - int j= ii; - while (!(*on)[j].Required()) { - Adj(on,j); // next geom edge - j=1-j; - if (e[Gh.Number(on)]) break; // optimisation - e[Gh.Number(on)] = ei; - } - } - - int kk=0; - for ( i=0;i2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield); + + if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) { + if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); + continue; + } + + double *sf = ss; + for (Int4 nufield=0;nufield2) printf(" Filling holes in mesh of %i vertices\n",nbv); - - //initialize ordre - if (!ordre){ - throw ErrorException(__FUNCT__,exprintf("!ordre")); - } - for (i=0;iaddtrie(Number(edges[i][0]),Number(edges[i][1]))); - } - if (kk != nbe) { - throw ErrorException(__FUNCT__,exprintf("Some Double edge in the mesh, the number is %i",kk-nbe)); - } - - // - for (i=0;iaddtrie(Number(triangles[i][VerticesOfTriangularEdge[j][0]]), - Number(triangles[i][VerticesOfTriangularEdge[j][1]])); - Int4 invisible = triangles[i].Hidden(j); - if(st[k]==-1){ - st[k]=3*i+j; - } - else if(st[k]>=0) { - if (triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3))){ - throw ErrorException(__FUNCT__,exprintf("(triangles[i].TriangleAdj(j) || triangles[st[k] / 3].TriangleAdj((int) (st[k]%3)))")); - } - - triangles[i].SetAdj2(j,triangles + st[k] / 3,(int) (st[k]%3)); - if (invisible) triangles[i].SetHidden(j); - if (k5) { - printf(" info on Mesh %s:\n",name); - printf(" - number of vertices = %i \n",nbv); - printf(" - number of triangles = %i \n",nbt); - printf(" - number of given edges = %i \n",nbe); - printf(" - number of all edges = %i \n" ,edge4->nb()); - printf(" - Euler number 1 - nb of holes = %i \n" ,nbt-edge4->nb()+nbv); - } - - // check the consistant of edge[].adj and the geometrical required vertex - Int4 k=0; - for (i=0;inb();i++){ - if (st[i] >=0){ // edge alone - if (i < nbe) { - Int4 i0=edge4->i(i); - ordre[i0] = vertices+i0; - Int4 i1=edge4->j(i); - ordre[i1] = vertices+i1; - } - else { - k=k+1; - if (k <20) { - throw ErrorException(__FUNCT__,exprintf("Lost boundary edges %i : %i %i",i,edge4->i(i),edge4->j(i))); - } - } - } - } - if(k != 0) { - throw ErrorException(__FUNCT__,exprintf("%i boundary edges are not defined as edges",k)); - } - - /* mesh generation with boundary points*/ - Int4 nbvb = 0; - for (i=0;ii, ordre[1]->i, ordre[i]->i ) == 0;) - if ( ++i >= nbvb) { - throw ErrorException(__FUNCT__,exprintf("FillHoleInMesh: All the vertices are aligned")); - } - Exchange( ordre[2], ordre[i]); - - Vertex * v0=ordre[0], *v1=ordre[1]; - - - triangles[0](0) = 0; // sommet pour infini - triangles[0](1) = v0; - triangles[0](2) = v1; - - triangles[1](0) = 0;// sommet pour infini - triangles[1](2) = v0; - triangles[1](1) = v1; - const int e0 = OppositeEdge[0]; - const int e1 = NextEdge[e0]; - const int e2 = PreviousEdge[e0]; - triangles[0].SetAdj2(e0, &triangles[1] ,e0); - triangles[0].SetAdj2(e1, &triangles[1] ,e2); - triangles[0].SetAdj2(e2, &triangles[1] ,e1); - - triangles[0].det = -1; // faux triangles - triangles[1].det = -1; // faux triangles - - triangles[0].SetTriangleContainingTheVertex(); - triangles[1].SetTriangleContainingTheVertex(); - - triangles[0].link=&triangles[1]; - triangles[1].link=&triangles[0]; - - if (!quadtree ) - delete quadtree; // ->ReInitialise(); - - quadtree = new QuadTree(this,0); - quadtree->Add(*v0); - quadtree->Add(*v1); - - // We add the vertices one by one - Int4 NbSwap=0; - for (Int4 icount=2; icounti,dete); - quadtree->Add(*vi); - Add(*vi,tcvi,dete); - NbSwap += vi->Optim(1,1); - } - - // inforce the boundary - TriangleAdjacent ta(0,0); - Int4 nbloss = 0,knbe=0; - for ( i = 0; i < nbe; i++){ - if (st[i] >=0){ // edge alone => on border ... FH oct 2009 - Vertex & a=edges[i][0], & b = edges[i][1]; - if (a.t && b.t) // le bug est la si maillage avec des bod non raffine 1. - { - knbe++; - if (ForceEdge(a,b,ta)<0) - nbloss++; - } - } - } - if(nbloss) { - throw ErrorException(__FUNCT__,exprintf("we lost(?) %i edges other %i",nbloss,knbe)); - } - - FindSubDomain(1); - // remove all the hole - // remove all the good sub domain - Int4 krm =0; - for (i=0;iaddtrie(v0?Number(v0):nbv,v1? Number(v1):nbv); - if (st[k]<0){ - throw ErrorException(__FUNCT__,exprintf("st[k]<0")); - } - tta.SetAdj2(ja,savetriangles + st[k] / 3,(int) (st[k]%3)); - ta.SetLock(); - st[k]=-2-st[k]; - } - } - } - } - Int4 NbTfillHoll =0; - for (i=0;isavenbtx ){ - throw ErrorException(__FUNCT__,exprintf("savenbt+NbTfillHoll>savenbtx")); - } - // copy of the outside triangles in saveTriangles - for (i=0;i=0) { - savetriangles[savenbt]=triangles[i]; - savetriangles[savenbt].link=0; - savenbt++; - } - } - // gestion of the adj - k =0; - Triangle * tmax = triangles + nbt; - for (i=0;i= triangles && ta < tmax) + int ie = OppositeEdge[j]; + TriangleAdjacent ta = triangles[i].Adj(ie); + Triangle *tt = ta; + if (tt && tt->link) { - ta= savetriangles + ta->color; - ti.SetAdj2(j,ta,aa); - if(lck) ti.SetLocked(j); + Vertex &v = *ta.OppositeVertex(); + R2 V = v; + Int4 iV = Number(v); + Real8 lA = bamg::Area2(V,B,C)/dd; + Real8 lB = bamg::Area2(A,V,C)/dd; + Real8 lC = bamg::Area2(A,B,V)/dd; + taa[0][j] = lB*lC; + taa[1][j] = lC*lA; + taa[2][j] = lA*lB; + lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ; + + bb[j] = ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ; + } + else + { + nbb++; + taa[0][j]=0; + taa[1][j]=0; + taa[2][j]=0; + taa[j][j]=1; + bb[j]=0; } } + + // resolution of 3x3 lineaire system transpose + Real8 det33 = det3x3(taa[0],taa[1],taa[2]); + Real8 cBC = det3x3(bb,taa[1],taa[2]); + Real8 cCA = det3x3(taa[0],bb,taa[2]); + Real8 cAB = det3x3(taa[0],taa[1],bb); + + if (!det33){ + throw ErrorException(__FUNCT__,exprintf("!det33")); + } + // computation of the gradient in the element + + // H( li*lj) = grad li grad lj + grad lj grad lj + // grad li = njk / detT ; with i j k =(A,B,C) + Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x); + Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y); + Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) + + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y); + Real8 coef = 1.0/(3*dd*det33); + Real8 coef2 = 2*coef; + Hxx *= coef2; + Hyy *= coef2; + Hxy *= coef2; + if(nbb==0) + { + dxdx[iA] += Hxx; + dydy[iA] += Hyy; + dxdy[iA] += Hxy; + + dxdx[iB] += Hxx; + dydy[iB] += Hyy; + dxdy[iB] += Hxy; + + dxdx[iC] += Hxx; + dydy[iC] += Hyy; + dxdy[iC] += Hxy; + } + } - // OutSidesTriangles = triangles; - // Int4 NbOutSidesTriangles = nbt; - - // restore triangles; - nbt=savenbt; - nbtx=savenbtx; - delete [] triangles; - delete [] subdomains; - triangles = savetriangles; - subdomains = savesubdomains; - if (k) { - throw ErrorException(__FUNCT__,exprintf("number of triangles edges alone = %i",k)); - } - FindSubDomain(); - - delete edge4; - delete [] st; - for (i=0;iAdd(vertices[i]); - - SetVertexFieldOn(); - - for (i=0;iIsRequiredVertex()) { - throw ErrorException(__FUNCT__,exprintf("adj and vertex required esge(?)")); - } + else + { + + // if edge on boundary no contribution => normal = 0 + if ( ! tBC || ! tBC->link ) nBC = O; + if ( ! tCA || ! tCA->link ) nCA = O; + if ( ! tAB || ! tAB->link ) nAB = O; + + // remark we forgot a 1/2 because + // $\\int_{edge} w_i = 1/2 $ if $i$ is in edge + // 0 if not + // if we don't take the boundary + // dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; + + dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; + dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x; + dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x; + + // warning optimization (1) the divide by 2 is done on the metrix construction + dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ; + dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ; + dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; + + dydy[iA] += ( nCA.y + nAB.y ) *Grads.y; + dydy[iB] += ( nAB.y + nBC.y ) *Grads.y; + dydy[iC] += ( nBC.y + nCA.y ) *Grads.y; + } + + } // for real all triangles + Int4 kk=0; + for ( iv=0,k=0 ; iv0) + { + dxdx[iv] /= 2*Mmassxx[iv]; + // warning optimization (1) on term dxdy[iv]*ci/2 + dxdy[iv] /= 4*Mmassxx[iv]; + dydy[iv] /= 2*Mmassxx[iv]; + // Compute the matrix with abs(eigen value) + Metric M(dxdx[iv], dxdy[iv], dydy[iv]); + MatVVP2x2 Vp(M); + Vp.Abs(); + M = Vp; + dxdx[iv] = M.a11; + dxdy[iv] = M.a21; + dydy[iv] = M.a22; + } + else kk++; + + + // correction of second derivate + // by a laplacien + + Real8 *d2[3] = { dxdx, dxdy, dydy}; + Real8 *dd; + for (int xy = 0;xy<3;xy++) + { + dd = d2[xy]; + // do leat 2 iteration for boundary problem + for (int ijacobi=0;ijacobi0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); + workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; + } + for (iv=0;iv2) { + printf(" Field %i of solution %i\n",nufield,nusol); + printf(" before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5)); + printf(" after bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); + } + } // end of for all field + }// end for all solution + + delete [] detT; + delete [] Mmass; + delete [] dxdx; + delete [] dxdy; + delete [] dydy; + delete [] workT; + delete [] workV; + delete [] Mmassxx; + delete [] OnBoundary; + +} +/*}}}1*/ +/*FUNCTION Triangles::IntersectGeomMetric{{{1*/ +void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0){ + long int verbosity=0; + Real8 ss[2]={0.00001,0.99999}; + Real8 errC = 2*sqrt(2*err); + Real8 hmax = Gh.MaximalHmax(); + Real8 hmin = Gh.MinimalHmin(); + Real8 maxaniso = 1e6; + if (hmax<=0){ + throw ErrorException(__FUNCT__,exprintf("hmax<=0")); + } + SetVertexFieldOn(); + if (errC > 1) errC = 1; + for (Int4 i=0;iR1tg(s,tg); + Real8 ht = hmax; + if (R1>1.0e-20) { // err relative to the length of the edge + ht = Min(Max(errC/R1,hmin),hmax); + } + Real8 hn = iso? ht : Min(hmax,ht*maxaniso); + if (ht<=0 || hn<=0){ + throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0")); + } + MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg); + Metric MVp(Vp); + edges[i][j].m.IntersectWith(MVp); + } + + } + // the problem is for the vertex on vertex +} +/*}}}1*/ + /*FUNCTION Triangles::Insert{{{1*/ + void Triangles::Insert() { + long int verbosity=2; + if (verbosity>2) printf(" Insert initial %i vertices\n",nbv); + Triangles * OldCurrentTh =CurrentTh; + + CurrentTh=this; + SetIntCoor(); + Int4 i; + for (i=0;ii, ordre[1]->i, ordre[i]->i ) == 0;){ + if ( ++i >= nbv) { + throw ErrorException(__FUNCT__,exprintf("all the vertices are aligned")); + } + } + // echange i et 2 dans ordre afin + // que les 3 premiers ne soit pas aligne + Exchange( ordre[2], ordre[i]); + + // on ajoute un point a l'infini pour construire le maillage + // afin d'avoir une definition simple des aretes frontieres + nbt = 2; + + // on construit un maillage trivale forme + // d'une arete et de 2 triangles + // construit avec le 2 aretes orientes et + Vertex * v0=ordre[0], *v1=ordre[1]; + + triangles[0](0) = 0; // sommet pour infini + triangles[0](1) = v0; + triangles[0](2) = v1; + + triangles[1](0) = 0;// sommet pour infini + triangles[1](2) = v0; + triangles[1](1) = v1; + const int e0 = OppositeEdge[0]; + const int e1 = NextEdge[e0]; + const int e2 = PreviousEdge[e0]; + triangles[0].SetAdj2(e0, &triangles[1] ,e0); + triangles[0].SetAdj2(e1, &triangles[1] ,e2); + triangles[0].SetAdj2(e2, &triangles[1] ,e1); + + triangles[0].det = -1; // faux triangles + triangles[1].det = -1; // faux triangles + + triangles[0].SetTriangleContainingTheVertex(); + triangles[1].SetTriangleContainingTheVertex(); + + triangles[0].link=&triangles[1]; + triangles[1].link=&triangles[0]; + + // nbtf = 2; + if ( !quadtree ) quadtree = new QuadTree(this,0); + quadtree->Add(*v0); + quadtree->Add(*v1); + + //the vertices are added one by one + Int4 NbSwap=0; + + + if (verbosity>3) printf(" Begining of insertion process...\n"); + + for (Int4 icount=2; icounti,dete); + quadtree->Add(*vi); + Add(*vi,tcvi,dete); + NbSwap += vi->Optim(1,0); + } + if (verbosity>3) { + printf(" NbSwap of insertion: %i\n",NbSwap); + printf(" NbSwap/nbv: %i\n",NbSwap/nbv); + printf(" NbUnSwap: %i\n",NbUnSwap); + printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); + } + NbUnSwap = 0; + // construction d'un ordre aleatoire + // const int PrimeNumber= (nbv % 999983) ? 1000003: 999983 ; +#ifdef NBLOOPOPTIM + + k3 = rand()%nbv ; + for (int is4=0; is4Optim(0,0); + if (verbosity>3) { + printf(" Optim Loop: %i\n",Nbloop); + printf(" NbSwap/nbv: %i\n",NbSwap/nbv); + printf(" NbUnSwap: %i\n",NbUnSwap); + printf(" NbUnSwap/nbv %i\n",NbUnSwap/nbv); + } + NbUnSwap = 0; + if(!NbSwap) break; } + ReMakeTriangleContainingTheVertex(); + // because we break the TriangleContainingTheVertex +#endif CurrentTh=OldCurrentTh; } /*}}}1*/ - /*FUNCTION Triangles::Optim{{{1*/ - Int4 Triangle::Optim(Int2 i,int koption) { - // turne around in positif sens - Int4 NbSwap =0; - Triangle *t = this; - int k=0,j =OppositeEdge[i]; - int jp = PreviousEdge[j]; - // initialise tp, jp the previous triangle & edge - Triangle *tp= at[jp]; - jp = aa[jp]&3; + /*FUNCTION Triangles::InsertNewPoints{{{1*/ + Int4 Triangles::InsertNewPoints(Int4 nbvold,Int4 & NbTSwap) { + long int verbosity=2; + Real8 seuil= 1.414/2 ;// for two close point + Int4 i; + // insertion part --- + + const Int4 nbvnew = nbv-nbvold; + if (verbosity>5) printf(" Try to Insert the %i new points\n",nbvnew); + Int4 NbSwap=0; + Icoor2 dete[3]; + + // construction d'un ordre aleatoire + if (! nbvnew) + return 0; + if (nbvnew) { + const Int4 PrimeNumber= AGoodNumberPrimeWith(nbv) ; + Int4 k3 = rand()%nbvnew ; + for (Int4 is3=0; is3ReferenceNumber=i; + } + // be carefull + Int4 iv = nbvold; + for (i=nbvold;iToClose(vi,seuil,hi,hj)) + { + // a good new point + Vertex & vj = vertices[iv]; + Int4 j = vj.ReferenceNumber; + if ( &vj!= ordre[j]){ + throw ErrorException(__FUNCT__,exprintf("&vj!= ordre[j]")); + } + if(i!=j) + { // for valgring + Exchange(vi,vj); + Exchange(ordre[j],ordre[i]); + } + vj.ReferenceNumber=0; + Triangle *tcvj = FindTriangleContening(vj.i,dete); + if (tcvj && !tcvj->link){ + throw ErrorException(__FUNCT__,exprintf("problem inserting point")); + } + quadtree->Add(vj); + Add(vj,tcvj,dete); + NbSwap += vj.Optim(1); + iv++; + } + } + if (verbosity>3) { + printf(" number of new points: %i\n",iv); + printf(" number of to close (?) points: %i\n",nbv-iv); + printf(" number of swap after: %i\n",NbSwap); + } + nbv = iv; + } + + for (i=nbvold;i3) printf(" NbSwap=%i\n",NbSwap); + + NbTSwap += NbSwap ; + return nbv-nbvold; + } + /*}}}1*/ +/*FUNCTION Triangles::MakeGeometricalEdgeToEdge{{{1*/ +Edge** Triangles::MakeGeometricalEdgeToEdge() { + if (!Gh.nbe){ + throw ErrorException(__FUNCT__,exprintf("!Gh.nbe")); + } + Edge **e= new (Edge* [Gh.nbe]); + + Int4 i; + for ( i=0;ion; + e[Gh.Number(on)] = ei; + } + for ( i=0;ion; + int j= ii; + while (!(*on)[j].Required()) { + Adj(on,j); // next geom edge + j=1-j; + if (e[Gh.Number(on)]) break; // optimisation + e[Gh.Number(on)] = ei; + } + } + + int kk=0; + for ( i=0;i2) printf("MakeQuadrangles costheta = %g\n",costheta); + + if (costheta >1) { + if (verbosity>5) printf(" do nothing: costheta > 1\n"); + } + + + Int4 nbqq = (nbt*3)/2; + DoubleAndInt4 *qq = new DoubleAndInt4[nbqq]; + + Int4 i,ij; + int j; + Int4 k=0; + for (i=0;i=costheta) + qq[k++].i3j=i*3+j; + // sort qq + HeapSort(qq,k); + + Int4 kk=0; + for (ij=0;ij=costheta) + triangles[i].SetHidden(j),kk++; + } + NbOfQuad = kk; + if (verbosity>2){ + printf(" number of quadrilaterals = %i\n",NbOfQuad); + printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); + printf(" number of outside triangles = %i\n",NbOutT); + } + delete [] qq; + } + /*}}}1*/ +/*FUNCTION Triangles::MakeQuadTree{{{1*/ +void Triangles::MakeQuadTree() { + long int verbosity=0; + if ( !quadtree ) quadtree = new QuadTree(this); + +} +/*}}}1*/ +/*FUNCTION Triangles::MaxSubDivision{{{1*/ +void Triangles::MaxSubDivision(Real8 maxsubdiv) { + long int verbosity=0; + + const Real8 maxsubdiv2 = maxsubdiv*maxsubdiv; + if(verbosity>1) printf(" Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv); + // for all the edges + // if the len of the edge is to long + Int4 it,nbchange=0; + Real8 lmax=0; + for (it=0;it maxsubdiv2) + { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M + Real8 lc = M(AC,AC); + D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; + D2xD2 Rt1(Rt.inv()); + D2xD2 D(maxsubdiv2,0,0,lc); + D2xD2 MM = Rt1*D*Rt1.t(); + v0.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y); + nbchange++; + } + M = v1; + l = M(AB,AB); + lmax = Max(lmax,l); + if(l> maxsubdiv2) + { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M + Real8 lc = M(AC,AC); + D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; + D2xD2 Rt1(Rt.inv()); + D2xD2 D(maxsubdiv2,0,0,lc); + D2xD2 MM = Rt1*D*Rt1.t(); + v1.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y); + nbchange++; + } + + + } + } + } + if(verbosity>3){ + printf(" number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5)); + } +} +/*}}}1*/ + /*FUNCTION Triangles::MetricAt{{{1*/ + Metric Triangles::MetricAt(const R2 & A) const { + I2 a = toI2(A); + Icoor2 deta[3]; + Triangle * t =FindTriangleContening(a,deta); + if (t->det <0) { // outside + double ba,bb; + TriangleAdjacent edge= CloseBoundaryEdge(a,t,ba,bb) ; + return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));} + else { // inside + Real8 aa[3]; + Real8 s = deta[0]+deta[1]+deta[2]; + aa[0]=deta[0]/s; + aa[1]=deta[1]/s; + aa[2]=deta[2]/s; + return Metric(aa,(*t)[0],(*t)[1],(*t)[2]); + } + } + /*}}}1*/ +/*FUNCTION Triangles::NearestVertex{{{1*/ +Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j) { + return quadtree->NearestVertex(i,j); +} +/*}}}1*/ + /*FUNCTION Triangles::NewPoints{{{1*/ + void Triangles::NewPoints(Triangles & Bh,int KeepBackVertex) { + long int verbosity=2; + Int4 nbtold(nbt),nbvold(nbv); + if (verbosity>2) printf(" Triangles::NewPoints\n"); + Int4 i,k; + int j; + Int4 *first_np_or_next_t = new Int4[nbtx]; + Int4 NbTSwap =0; + // insert old point + nbtold = nbt; + + if (KeepBackVertex && (&Bh != this) && (nbv+Bh.nbv< nbvx)){ + // Bh.SetVertexFieldOn(); + for (i=0;i= nbt + // the list of test triangle is + // the next traingle on i is -first_np_or_next_t[i] + int iter=0; + // Big loop do { - while (t->swap(j,koption)) - { - NbSwap++; - k++; - if (k>=20000){ - throw ErrorException(__FUNCT__,exprintf("k>=20000")); + iter++; + nbtold = nbt; + nbvold = nbv; + + // default size of IntersectionTriangle + + i=Headt; + next_t=-first_np_or_next_t[i]; + for(t=&triangles[i];i=nbt ){ + throw ErrorException(__FUNCT__,exprintf("i<0 || i>=nbt")); } - t= tp->at[jp]; // set unchange t qnd j for previous triangles - j= NextEdge[tp->aa[jp]&3]; - } - // end on this Triangle - tp = t; - jp = NextEdge[j]; - + first_np_or_next_t[i] = iter; + for(j=0;j<3;j++) + { // for each edge + TriangleAdjacent tj(t,j); + Vertex & vA = * tj.EdgeVertex(0); + Vertex & vB = * tj.EdgeVertex(1); + + if (!t->link) continue;// boundary + if (t->det <0) continue; + if (t->Locked(j)) continue; + + TriangleAdjacent tadjj = t->Adj(j); + Triangle * ta= tadjj; + + if (ta->det <0) continue; + + R2 A = vA; + R2 B = vB; + + k=Number(ta); + + if(first_np_or_next_t[k]==iter) // this edge is done before + continue; // next edge of the triangle + + //const Int4 NbvOld = nbv; + lIntTria.SplitEdge(Bh,A,B); + lIntTria.NewPoints(vertices,nbv,nbvx); + } // end loop for each edge + + }// for triangle + + if (!InsertNewPoints(nbvold,NbTSwap)) + break; + + for (i=nbtold;it, EdgesVertexTriangle[s->vint][1]); + Triangle * tbegin= (Triangle*) ta; + Int4 kt; + do { + kt = Number((Triangle*) ta); + if (first_np_or_next_t[kt]>0) + first_np_or_next_t[kt]=-Headt,Headt=kt; + if (ta.EdgeVertex(0)!=s){ + throw ErrorException(__FUNCT__,exprintf("ta.EdgeVertex(0)!=s")); + } + ta = Next(Adj(ta)); + } while ( (tbegin != (Triangle*) ta)); + } + + } while (nbv!=nbvold); + + delete [] first_np_or_next_t; + + Int4 NbSwapf =0,NbSwp; + + // bofbof + + + NbSwp = NbSwapf; + for (i=0;iswap(j,koption)) + { + NbSwap++; + k++; + if (k>=20000){ + throw ErrorException(__FUNCT__,exprintf("k>=20000")); + } t= tp->at[jp]; // set unchange t qnd j for previous triangles j= NextEdge[tp->aa[jp]&3]; - - } while( t != this); - return NbSwap; + } + // end on this Triangle + tp = t; + jp = NextEdge[j]; + + t= tp->at[jp]; // set unchange t qnd j for previous triangles + j= NextEdge[tp->aa[jp]&3]; + + } while( t != this); + return NbSwap; +} +/*}}}1*/ +/*FUNCTION Triangles::PreInit{{{1*/ +void Triangles::PreInit(Int4 inbvx,char *fname) { + long int verbosity=0; + + srand(19999999); + NbRef=0; + // allocGeometry=0; + identity=0; + NbOfTriangleSearchFind =0; + NbOfSwapTriangle =0; + nbiv=0; + nbv=0; + nbvx=inbvx; + nbt=0; + NbOfQuad = 0; + nbtx=2*inbvx-2; + NbSubDomains=0; + NbVertexOnBThVertex=0; + NbVertexOnBThEdge=0; + VertexOnBThVertex=0; + VertexOnBThEdge=0; + + NbCrackedVertices=0; + NbCrackedEdges =0; + CrackedEdges =0; + nbe = 0; + name = fname ; + + if (inbvx) { + vertices=new Vertex[nbvx]; + if (!vertices){ + throw ErrorException(__FUNCT__,exprintf("!vertices")); + } + ordre=new (Vertex* [nbvx]); + if (!ordre){ + throw ErrorException(__FUNCT__,exprintf("!ordre")); + } + triangles=new Triangle[nbtx]; + if (!triangles){ + throw ErrorException(__FUNCT__,exprintf("!triangles")); + } + } + else { + vertices=0; + ordre=0; + triangles=0; + nbtx=0; + } + if ( name || inbvx) { + time_t timer =time(0); + char buf[70]; + strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer)); + counter++; + char countbuf[30]; + sprintf(countbuf,"%d",counter); + int lg =0 ; + if (&BTh != this && BTh.name) + lg = strlen(BTh.name)+4; + identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2 + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)]; + identity[0]=0; + if (lg) + strcat(strcat(strcat(identity,"B="),BTh.name),", "); + + if (Gh.name) + strcat(strcat(identity,"G="),Gh.name); + strcat(strcat(identity,";"),countbuf); + strcat(identity,buf); + } + + quadtree=0; + edges=0; + VerticesOnGeomVertex=0; + VerticesOnGeomEdge=0; + NbVerticesOnGeomVertex=0; + NbVerticesOnGeomEdge=0; + subdomains=0; + NbSubDomains=0; +} +/*}}}1*/ + /*FUNCTION Triangles::ProjectOnCurve{{{1*/ + GeometricalEdge* Triangles::ProjectOnCurve( Edge & BhAB, Vertex & vA, Vertex & vB, + Real8 theta,Vertex & R,VertexOnEdge & BR,VertexOnGeom & GR) { + void *pA=0,*pB=0; + Real8 tA=0,tB=0; + R2 A=vA,B=vB; + Vertex * pvA=&vA, * pvB=&vB; + if (vA.vint == IsVertexOnVertex){ + pA=vA.onbv; + } + else if (vA.vint == IsVertexOnEdge){ + pA=vA.onbe->be; + tA=vA.onbe->abcisse; + } + else { + throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve On Vertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vA))); + } + + if (vB.vint == IsVertexOnVertex){ + pB=vB.onbv; + } + else if(vB.vint == IsVertexOnEdge){ + pB=vB.onbe->be; + tB=vB.onbe->abcisse; + } + else { + throw ErrorException(__FUNCT__,exprintf("ProjectOnCurve On Vertex %i forget call to SetVertexFieldOnBTh",BTh.Number(vB))); + } + Edge * e = &BhAB; + if (!pA || !pB || !e){ + throw ErrorException(__FUNCT__,exprintf("!pA || !pB || !e")); + } + // be carefull the back ground edge e is on same geom edge + // of the initiale edge def by the 2 vertex A B; + //check Is a background Mesh; + if (e=BTh.edges+BTh.nbe){ + throw ErrorException(__FUNCT__,exprintf("e=BTh.edges+BTh.nbe")); + } + // walk on BTh edge + //not finish ProjectOnCurve with BackGround Mesh); + // 1 first find a back ground edge contening the vertex A + // 2 walk n back gound boundary to find the final vertex B + + if( vA.vint == IsVertexOnEdge) + { // find the start edge + e = vA.onbe->be; + + } + else if (vB.vint == IsVertexOnEdge) + { + theta = 1-theta; + Exchange(tA,tB); + Exchange(pA,pB); + Exchange(pvA,pvB); + Exchange(A,B); + e = vB.onbe->be; + + } + else{ // do the search by walking + throw ErrorException(__FUNCT__,exprintf("case not supported yet")); + } + + // find the direction of walking with sens of edge and pA,PB; + R2 AB=B-A; + + Real8 cosE01AB = (( (R2) (*e)[1] - (R2) (*e)[0] ) , AB); + int kkk=0; + int sens = (cosE01AB>0) ? 1 : 0; + + // Real8 l=0; // length of the edge AB + Real8 abscisse = -1; + + for (int cas=0;cas<2;cas++){ + // 2 times algo: + // 1 for computing the length l + // 2 for find the vertex + int iii; + Vertex *v0=pvA,*v1; + Edge *neee,*eee; + Real8 lg =0; // length of the curve + Real8 te0; + // we suppose take the curve's abcisse + for ( eee=e,iii=sens,te0=tA; + eee && ((( void*) eee) != pB) && (( void*) (v1=&((*eee)[iii]))) != pB ; + neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) { + + kkk=kkk+1; + if (kkk>=100){ + throw ErrorException(__FUNCT__,exprintf("kkk>=100")); + } + if (!eee){ + throw ErrorException(__FUNCT__,exprintf("!eee")); + } + Real8 lg0 = lg; + Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); + lg += dp; + if (cas && abscisse <= lg) { // ok we find the geom edge + Real8 sss = (abscisse-lg0)/dp; + Real8 thetab = te0*(1-sss)+ sss*iii; + if (thetab<0 || thetab>1){ + throw ErrorException(__FUNCT__,exprintf("thetab<0 || thetab>1")); + } + BR = VertexOnEdge(&R,eee,thetab); + return Gh.ProjectOnCurve(*eee,thetab,R,GR); + } + } + // we find the end + if (v1 != pvB){ + if (( void*) v1 == pB) + tB = iii; + + Real8 lg0 = lg; + if (!eee){ + throw ErrorException(__FUNCT__,exprintf("!eee")); + } + v1 = pvB; + Real8 dp = LengthInterpole(v0->m,v1->m,(R2) *v1 - (R2) *v0); + lg += dp; + abscisse = lg*theta; + if (abscisse <= lg && abscisse >= lg0 ) // small optimisation we know the lenght because end + { // ok we find the geom edge + Real8 sss = (abscisse-lg0)/dp; + Real8 thetab = te0*(1-sss)+ sss*tB; + if (thetab<0 || thetab>1){ + throw ErrorException(__FUNCT__,exprintf("thetab<0 || thetab>1")); + } + BR = VertexOnEdge(&R,eee,thetab); + return Gh.ProjectOnCurve(*eee,thetab,R,GR); + } + } + abscisse = lg*theta; + + } + throw ErrorException(__FUNCT__,exprintf("Big bug...")); + return 0; // just for the compiler + } + /*}}}1*/ + /*FUNCTION Triangles::ReNumberingTheTriangleBySubDomain{{{1*/ + void Triangles::ReNumberingTheTriangleBySubDomain(bool justcompress) { + long int verbosity=0; + Int4 *renu= new Int4[nbt]; + register Triangle *t0,*t,*te=triangles+nbt; + register Int4 k=0,it,i,j; + + for ( it=0;it= nbt ){ + throw ErrorException(__FUNCT__,exprintf("kt<0 || kt >= nbt")); + } + if (renu[kt]!=-1){ + throw ErrorException(__FUNCT__,exprintf("renu[kt]!=-1")); + } + renu[kt]=k++; + } + while (t0 != (t=t->link)); + } + // take is same numbering if possible + if(justcompress) + for ( k=0,it=0;it=0 ) + renu[it]=k++; + + // put the outside triangles at the end + for ( it=0;it= 0) // a new sub cycle + { + i=it; + Triangle ti=triangles[i],tj; + while ( (j=renu[i]) >= 0) + { // i is old, and j is new + renu[i] = -1; // mark + tj = triangles[j]; // save new + triangles[j]= ti; // new <- old + i=j; // next + ti = tj; + } + } + delete [] renu; + nt = nbt - NbOutT; + } /*}}}1*/ - /*FUNCTION Triangles::SmoothingVertex{{{1*/ - void Triangles::SmoothingVertex(int nbiter,Real8 omega ) { - long int verbosity=0; - // if quatree exist remove it end reconstruct - if (quadtree) delete quadtree; - quadtree=0; - ReMakeTriangleContainingTheVertex(); - Triangle vide; // a triangle to mark the boundary vertex - Triangle ** tstart= new Triangle* [nbv]; - Int4 i,j,k; - // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide - if ( this == & BTh) - for ( i=0;i2) printf(" SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega); - for (k=0;k3) printf(" move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap); + Vertex *v = VerticesOnGeomVertex[i].mv; + if (v>=vertices && v < ve) + VerticesOnGeomVertex[i].mv=vertices+renu[Number(v)]; } - delete [] tstart; - if (quadtree) quadtree= new QuadTree(this); + for (i=0;i< NbVerticesOnGeomEdge;i++) + { + Vertex *v =VerticesOnGeomEdge[i].mv; + if (v>=vertices && v < ve) + VerticesOnGeomEdge[i].mv=vertices+renu[Number(v)]; + } + + for (i=0;i< NbVertexOnBThVertex;i++) + { + Vertex *v=VertexOnBThVertex[i].v; + if (v>=vertices && v < ve) + VertexOnBThVertex[i].v=vertices+renu[Number(v)]; + } + + for (i=0;i< NbVertexOnBThEdge;i++) + { + Vertex *v=VertexOnBThEdge[i].v; + if (v>=vertices && v < ve) + VertexOnBThEdge[i].v=vertices+renu[Number(v)]; + } + + // move the Vertices without a copy of the array + // be carefull not trivial code + Int4 j; + for ( it=0;it= 0) // a new sub cycle + { + i=it; + Vertex ti=vertices[i],tj; + while ( (j=renu[i]) >= 0) + { // i is old, and j is new + renu[i] = -1-renu[i]; // mark + tj = vertices[j]; // save new + vertices[j]= ti; // new <- old + i=j; // next + ti = tj; + } + } + if (quadtree) + { delete quadtree; + quadtree = new QuadTree(this); + } + for ( it=0;itlink || Number(ta) >= it) + { + Vertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]]; + Vertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]]; + if ( !& vP || !&vQ) continue; + R2 PQ = vQ.r - vP.r; + Real8 l = log(LengthInterpole(vP,vQ,PQ)); + nbedges++; + k = (int) ((l - lmin)*delta); + k = Min(Max(k,0L),kmax); + histo[k]++; + } + } + } + printf(" --- Histogram of the unit mesh, nb of edges = %i\n",nbedges); + printf(" length of edge in | %% of edge | Nb of edges \n"); + printf(" --------------------+-------------+-------------\n"); + for (i=0;i<=kmax;i++){ + if (i==0) printf(" %10i",0); + else printf(" %10g",exp(lmin+i/delta)); + if (i==kmax) printf(" +inf "); + else printf(" %10g",exp(lmin+(i+1)/delta)); + printf("| %10g |\n",((long) ((10000.0 * histo[i])/ nbedges))/100.0); + printf(" %g\n",histo[i]); + } + printf(" --------------------+-------------+-------------\n"); +} +/*}}}1*/ +/*FUNCTION Triangles::SmoothingVertex{{{1*/ +void Triangles::SmoothingVertex(int nbiter,Real8 omega ) { + long int verbosity=0; + // if quatree exist remove it end reconstruct + if (quadtree) delete quadtree; + quadtree=0; + ReMakeTriangleContainingTheVertex(); + Triangle vide; // a triangle to mark the boundary vertex + Triangle ** tstart= new Triangle* [nbv]; + Int4 i,j,k; + // attention si Background == Triangle alors on ne peut pas utiliser la rechech rapide + if ( this == & BTh) + for ( i=0;i2) printf(" SmoothingVertex: nb Iteration = %i, Omega=%g\n",nbiter,omega); + for (k=0;k3) printf(" move max = %g, iteration = %i, nb of swap = %i\n",pow(delta,0.5),k,NbSwap); + } + + delete [] tstart; + if (quadtree) quadtree= new QuadTree(this); +} +/*}}}1*/ +/*FUNCTION Triangles::SmoothMetric{{{1*/ +void Triangles::SmoothMetric(Real8 raisonmax) { + long int verbosity=0; + + if(raisonmax<1.1) return; + if(verbosity > 1) printf(" Triangles::SmoothMetric raisonmax = %g\n",raisonmax); + ReMakeTriangleContainingTheVertex(); + Int4 i,j,kch,kk,ip; + Int4 *first_np_or_next_t0 = new Int4[nbv]; + Int4 *first_np_or_next_t1 = new Int4[nbv]; + Int4 Head0 =0,Head1=-1; + Real8 logseuil= log(raisonmax); + + for(i=0;i=0&& kk++<100) { + kch=0; + for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) { + // pour tous les triangles autour du sommet s + register Triangle* t= vertices[i].t; + if (!t){ + throw ErrorException(__FUNCT__,exprintf("!t")); + } + Vertex & vi = vertices[i]; + TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]); + Vertex *pvj0 = ta.EdgeVertex(0); + while (1) { + ta=Previous(Adj(ta)); + if (vertices+i != ta.EdgeVertex(1)){ + throw ErrorException(__FUNCT__,exprintf("vertices+i != ta.EdgeVertex(1)")); + } + Vertex & vj = *(ta.EdgeVertex(0)); + if ( &vj ) { + j= &vj-vertices; + if (j<0 || j >= nbv){ + throw ErrorException(__FUNCT__,exprintf("j<0 || j >= nbv")); + } + R2 Aij = (R2) vj - (R2) vi; + Real8 ll = Norme2(Aij); + if (0) { + Real8 hi = ll/vi.m(Aij); + Real8 hj = ll/vj.m(Aij); + if(hi < hj) + { + Real8 dh=(hj-hi)/ll; + if (dh>logseuil) { + vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi)); + if(first_np_or_next_t1[j]<0) + kch++,first_np_or_next_t1[j]=Head1,Head1=j; + } + } + } + else + { + Real8 li = vi.m(Aij); + if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) + if(first_np_or_next_t1[j]<0) // if the metrix change + kch++,first_np_or_next_t1[j]=Head1,Head1=j; + } + } + if ( &vj == pvj0 ) break; + } + } + Head0 = Head1; + Head1 = -1; + Exchange(first_np_or_next_t0,first_np_or_next_t1); + } + if(verbosity>2) printf(" number of iterations = %i\n",kch); + delete [] first_np_or_next_t0; + delete [] first_np_or_next_t1; +} +/*}}}1*/ + /*FUNCTION Triangles::SplitElement{{{1*/ + int Triangles::SplitElement(int choice){ long int verbosity=0; - if ( !quadtree ) quadtree = new QuadTree(this); - + + Direction NoDirOfSearch; + const int withBackground = &BTh != this && &BTh; + + ReNumberingTheTriangleBySubDomain(); + //int nswap =0; + const Int4 nfortria( choice ? 4 : 6); + if(withBackground) + { + BTh.SetVertexFieldOn(); + SetVertexFieldOnBTh(); + } + else + BTh.SetVertexFieldOn(); + + Int4 newnbt=0,newnbv=0; + Int4 * kedge = 0; + Int4 newNbOfQuad=NbOfQuad; + Int4 * ksplit= 0, * ksplitarray=0; + Int4 kkk=0; + int ret =0; + if (nbvxaddtrie(Number(edges[i][0]),Number(edges[i][1]))); + if (ong) // a geometrical edges + { + if (withBackground){ + // walk on back ground mesh + // newVertexOnBThEdge[ibe++] = VertexOnEdge(vertices[k],bedge,absicsseonBedge); + // a faire -- difficile + // the first PB is to now a background edge between the 2 vertices + if (!edgesGtoB){ + throw ErrorException(__FUNCT__,exprintf("!edgesGtoB")); + } + ong= ProjectOnCurve(*edgesGtoB[Gh.Number(edges[i].on)], + edges[i][0],edges[i][1],0.5,vertices[k], + newVertexOnBThEdge[kvb], + newVerticesOnGeomEdge[kvg++]); + vertices[k].ReferenceNumber= edges[i].ref; + vertices[k].DirOfSearch = NoDirOfSearch; + ; + // get the Info on background mesh + Real8 s = newVertexOnBThEdge[kvb]; + Vertex & bv0 = newVertexOnBThEdge[kvb][0]; + Vertex & bv1 = newVertexOnBThEdge[kvb][1]; + // compute the metrix of the new points + vertices[k].m = Metric(1-s,bv0,s,bv1); + kvb++; + } + else + { + ong=Gh.ProjectOnCurve(edges[i], + 0.5,vertices[k],newVerticesOnGeomEdge[kvg++]); + // vertices[k].i = toI2( vertices[k].r); + vertices[k].ReferenceNumber = edges[i].ref; + vertices[k].DirOfSearch = NoDirOfSearch; + vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); + } + } + else // straigth line edge --- + { + vertices[k].r = ((R2) edges[i][0] + (R2) edges[i][1] )*0.5; + vertices[k].m = Metric(0.5,edges[i][0],0.5,edges[i][1]); + vertices[k].on = 0; + } + //vertices[k].i = toI2( vertices[k].r); + R2 AB = vertices[k].r; + R2 AA = (A+AB)*0.5; + R2 BB = (AB+B)*0.5; + vertices[k].ReferenceNumber = edges[i].ref; + vertices[k].DirOfSearch = NoDirOfSearch; + + newedges[ie].on = Gh.Contening(AA,ong); + newedges[ie++].v[1]=vertices+k; + + newedges[ie]=edges[i]; + newedges[ie].adj[0]=newedges + ie -1; + newedges[ie].adj[1]=newedges+(edges[i].adj[1]-edges) ; + newedges[ie].on = Gh.Contening(BB,ong); + newedges[ie++].v[0]=vertices+k; + k++; + } + if (edgesGtoB) delete [] edgesGtoB; + edgesGtoB=0; + + newnbv=k; + newNbVerticesOnGeomEdge=kvg; + if (newnbv> nbvx) goto Error;// bug + + nbv = k; + + + kedge = new Int4[3*nbt+1]; + ksplitarray = new Int4[nbt+1]; + ksplit = ksplitarray +1; // because ksplit[-1] == ksplitarray[0] + + for (i=0;i<3*nbt;i++) + kedge[i]=-1; + + // + + for (i=0;i= lastT) + t.SetAdj2(j,0,0);// unset adj + const Vertex & v0 = t[VerticesOfTriangularEdge[j][0]]; + const Vertex & v1 = t[VerticesOfTriangularEdge[j][1]]; + Int4 ke =edge4->findtrie(Number(v0),Number(v1)); + if (ke>0) + { + Int4 ii = Number(tt); + int jj = ta; + Int4 ks = ke + nbvold; + kedge[3*i+j] = ks; + if (ii=0){ + // warning PB roundoff error + if (t.link && ( (aa=Area2( A.r , t[1].r , t[2].r )) < 0.0 + || (bb=Area2( t[0].r , A.r , t[2].r )) < 0.0 + || (cc=Area2( t[0].r , t[1].r , A.r )) < 0.0)){ + printf("%i not in triangle %i In= %i %i %i %i %i\n",ke + nbvold,i,!!t.link,aa,bb,cc,dd); + throw ErrorException(__FUNCT__,exprintf("Number of triangles with P2 interpolation Problem")); + } + } + else { + if (tt.link && ( (aa=Area2( A.r , tt[1].r , tt[2].r )) < 0 + || (bb=Area2( tt[0].r , A.r , tt[2].r )) < 0 + || (cc=Area2( tt[0].r , tt[1].r , A.r )) < 0)){ + printf("%i not in triangle %i In= %i %i %i %i %i\n",ke + nbvold,ii,!!tt.link,aa,bb,cc,dd); + throw ErrorException(__FUNCT__,exprintf("Number of triangles with P2 interpolation Problem")); + } + } + } + } + } + + for (i=0;ifindtrie(Number(v0),Number(v1)); + if (ke<0) // new + { + if (&tt) // internal triangles all the boundary + { // new internal edges + Int4 ii = Number(tt); + int jj = ta; + + kedge[3*i+j]=k;// save the vertex number + kedge[3*ii+jj]=k; + if (k=0 + kedge[3*i+j]=nbvold+ke; + kkk[nbsplitedge++]=j;// previously splited + } + } + else + kkk[nbsplitedge++]=j;// previously splited + + } + if (nbinvisible>=2){ + throw ErrorException(__FUNCT__,exprintf("nbinvisible>=2")); + } + switch (nbsplitedge) { + case 0: ksplit[i]=10; newnbt++; break; // nosplit + case 1: ksplit[i]=20+kkk[0];newnbt += 2; break; // split in 2 + case 2: ksplit[i]=30+3-kkk[0]-kkk[1];newnbt += 3; break; // split in 3 + case 3: + if (nbinvisible) ksplit[i]=40+invisibleedge,newnbt += 4; + else ksplit[i]=10*nfortria,newnbt+=nfortria; + break; + } + if (ksplit[i]<40){ + throw ErrorException(__FUNCT__,exprintf("ksplit[i]<40")); + } + } + // now do the element split + newNbOfQuad = 4*NbOfQuad; + nbv = k; + kkk = nbt; + ksplit[-1] = nbt; + // look on old true triangles + + for (i=0;i=7 || kk<=0){ + throw ErrorException(__FUNCT__,exprintf("kk>=7 || kk<=0")); + } + + // def the numbering k (edge) i vertex + int k0 = ke; + int k1 = NextEdge[k0]; + int k2 = PreviousEdge[k0]; + int i0 = OppositeVertex[k0]; + int i1 = OppositeVertex[k1]; + int i2 = OppositeVertex[k2]; + + Triangle &t0=triangles[i]; + Vertex * v0=t0(i0); + Vertex * v1=t0(i1); + Vertex * v2=t0(i2); + + if (nbmkadj>=10){ + throw ErrorException(__FUNCT__,exprintf("nbmkadj>=10")); + } + // -------------------------- + TriangleAdjacent ta0(t0.Adj(i0)),ta1(t0.Adj(i1)),ta2(t0.Adj(i2)); + // save the flag Hidden + int hid[]={t0.Hidden(0),t0.Hidden(1),t0.Hidden(2)}; + // un set all adj -- save Hidden flag -- + t0.SetAdj2(0,0,hid[0]); + t0.SetAdj2(1,0,hid[1]); + t0.SetAdj2(2,0,hid[2]); + // -- remake + switch (kk) { + case 1: break;// nothing + case 2: // + { + Triangle &t1=triangles[kkk++]; + t1=t0; + if (kedge[3*i+i0]<0){ + throw ErrorException(__FUNCT__,exprintf("kedge[3*i+i0]<0")); + } + Vertex * v3 = vertices + kedge[3*i+k0]; + + t0(i2) = v3; + t1(i1) = v3; + t0.SetAllFlag(k2,0); + t1.SetAllFlag(k1,0); + } + break; + case 3: // + { + Triangle &t1=triangles[kkk++]; + Triangle &t2=triangles[kkk++]; + t2=t1=t0; + if (kedge[3*i+k1]<0){ + throw ErrorException(__FUNCT__,exprintf("kedge[3*i+k1]<0")); + } + if (kedge[3*i+k2]<0){ + throw ErrorException(__FUNCT__,exprintf("kedge[3*i+k2]<0")); + } + + Vertex * v01 = vertices + kedge[3*i+k2]; + Vertex * v02 = vertices + kedge[3*i+k1]; + t0(i1) = v01; + t0(i2) = v02; + t1(i2) = v02; + t1(i0) = v01; + t2(i0) = v02; + t0.SetAllFlag(k0,0); + t1.SetAllFlag(k1,0); + t1.SetAllFlag(k0,0); + t2.SetAllFlag(k2,0); + } + break; + case 4: // + case 6: // split in 4 + { + Triangle &t1=triangles[kkk++]; + Triangle &t2=triangles[kkk++]; + Triangle &t3=triangles[kkk++]; + t3=t2=t1=t0; + if (kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0){ + throw ErrorException(__FUNCT__,exprintf("kedge[3*i+k0] <0 || kedge[3*i+k1]<0 || kedge[3*i+k2]<0")); + } + Vertex * v12 = vertices + kedge[3*i+k0]; + Vertex * v02 = vertices + kedge[3*i+k1]; + Vertex * v01 = vertices + kedge[3*i+k2]; + t0(i1) = v01; + t0(i2) = v02; + t0.SetAllFlag(k0,hid[k0]); + + t1(i0) = v01; + t1(i2) = v12; + t0.SetAllFlag(k1,hid[k1]); + + t2(i0) = v02; + t2(i1) = v12; + t2.SetAllFlag(k2,hid[k2]); + + t3(i0) = v12; + t3(i1) = v02; + t3(i2) = v01; + + t3.SetAllFlag(0,hid[0]); + t3.SetAllFlag(1,hid[1]); + t3.SetAllFlag(2,hid[2]); + + if ( kk == 6) + { + + Triangle &t4=triangles[kkk++]; + Triangle &t5=triangles[kkk++]; + + t4 = t3; + t5 = t3; + + t0.SetHidden(k0); + t1.SetHidden(k1); + t2.SetHidden(k2); + t3.SetHidden(0); + t4.SetHidden(1); + t5.SetHidden(2); + + if (nbv < nbvx ) + { + vertices[nbv].r = ((R2) *v01 + (R2) *v12 + (R2) *v02 ) / 3.0; + vertices[nbv].ReferenceNumber =0; + vertices[nbv].DirOfSearch =NoDirOfSearch; + //vertices[nbv].i = toI2(vertices[nbv].r); + Real8 a3[]={1./3.,1./3.,1./3.}; + vertices[nbv].m = Metric(a3,v0->m,v1->m,v2->m); + Vertex * vc = vertices +nbv++; + t3(i0) = vc; + t4(i1) = vc; + t5(i2) = vc; + + } + else + goto Error; + } + + } + break; + } + + // save all the new triangles + mkadj[nbmkadj++]=i; + Int4 jj; + if (t0.link) + for (jj=nbt;jj13){// 13 = 6 + 4 + + throw ErrorException(__FUNCT__,exprintf("nbmkadj>13")); + } + + if (kk==6) newNbOfQuad+=3; + for (jj=ksplit[i-1];jj= nbtx) goto Error; // bug + if (nbv>= nbvx) goto Error; // bug + // generation of the new triangles + + SetIntCoor("In SplitElement"); + + ReMakeTriangleContainingTheVertex(); + if(withBackground) + BTh.ReMakeTriangleContainingTheVertex(); + + delete [] edges; + edges = newedges; + nbe = newnbe; + NbOfQuad = newNbOfQuad; + + for (i=0;i2){ + printf(" number of quadrilaterals = %i\n",NbOfQuad); + printf(" number of triangles = %i\n",nbt-NbOutT- NbOfQuad*2); + printf(" number of outside triangles = %i\n",NbOutT); + } + + CurrentTh=OCurrentTh; + return 0; //ok + +Error: + nbv = nbvold; + nbt = nbtold; + NbOutT = NbOutTold; + // cleaning memory --- + delete newedges; + if (ksplitarray) delete [] ksplitarray; + if (kedge) delete [] kedge; + if (newVerticesOnGeomEdge) delete [] newVerticesOnGeomEdge; + if (edge4) delete edge4; + if(newVertexOnBThEdge) delete [] newVertexOnBThEdge; + + + CurrentTh= OCurrentTh; + return ret; // ok } /*}}}1*/ - /*FUNCTION Triangles::ShowRegulaty{{{1*/ - void Triangles::ShowRegulaty() const { - const Real8 sqrt32=sqrt(3.)*0.5; - const Real8 aireKh=sqrt32*0.5; - D2 Beq(1,0),Heq(0.5,sqrt32); - D2xD2 Br(D2xD2(Beq,Heq).t()); - D2xD2 B1r(Br.inv()); - double gammamn=1e100,hmin=1e100; - double gammamx=0,hmax=0; - double beta=1e100; - double beta0=0; - double alpha2=0; - double area=0,Marea=0; - // Real8 cf= Real8(coefIcoor); - // Real8 cf2= 6.*cf*cf; - int nt=0; - for (int it=0;itlink) { + printf("problem inserting point\n"); + } + + quadtree->Add(vi); + if (!tcvi || tcvi->det<0){// internal + throw ErrorException(__FUNCT__,exprintf("!tcvi || tcvi->det < 0")); + } + Add(vi,tcvi,dete); + NbSwap += vi.Optim(1); + iv++; + // } + } + if (verbosity>3) { + printf(" number of points: %i\n",iv); + printf(" number of swap to split internal edges with border vertices: %i\n",NbSwap); + nbv = iv; + } +} +if (NbSplitEdge>nbv-nbvold) printf("WARNING: not enough vertices to split all internal edges, we lost %i edges...\n",NbSplitEdge - ( nbv-nbvold)); +if (verbosity>2) printf("SplitInternalEdgeWithBorderVertices: Number of splited edge %i\n",NbSplitEdge); + +return NbSplitEdge; +} +/*}}}1*/ + /*FUNCTION Triangles::swap{{{1*/ + int Triangle::swap(Int2 a,int koption){ + if(a/4 !=0) return 0;// arete lock or MarkUnSwap + + register Triangle *t1=this,*t2=at[a];// les 2 triangles adjacent + register Int1 a1=a,a2=aa[a];// les 2 numero de l arete dans les 2 triangles + if(a2/4 !=0) return 0; // arete lock or MarkUnSwap + + register Vertex *sa=t1->ns[VerticesOfTriangularEdge[a1][0]]; + register Vertex *sb=t1->ns[VerticesOfTriangularEdge[a1][1]]; + register Vertex *s1=t1->ns[OppositeVertex[a1]]; + register Vertex *s2=t2->ns[OppositeVertex[a2]]; + + Icoor2 det1=t1->det , det2=t2->det ; + Icoor2 detT = det1+det2; + Icoor2 detA = Abs(det1) + Abs(det2); + Icoor2 detMin = Min(det1,det2); + + int OnSwap = 0; + // si 2 triangle infini (bord) => detT = -2; + if (sa == 0) {// les deux triangles sont frontieres + det2=bamg::det(s2->i,sb->i,s1->i); + OnSwap = det2 >0;} + else if (sb == 0) { // les deux triangles sont frontieres + det1=bamg::det(s1->i,sa->i,s2->i); + OnSwap = det1 >0;} + else if(( s1 != 0) && (s2 != 0) ) { + det1 = bamg::det(s1->i,sa->i,s2->i); + det2 = detT - det1; + OnSwap = (Abs(det1) + Abs(det2)) < detA; + + Icoor2 detMinNew=Min(det1,det2); + // if (detMin<0 && (Abs(det1) + Abs(det2) == detA)) OnSwap=BinaryRand();// just for test + if (! OnSwap &&(detMinNew>0)) { + OnSwap = detMin ==0; + if (! OnSwap) { + int kopt = koption; + while (1) + if(kopt) { + // critere de Delaunay pure isotrope + register Icoor2 xb1 = sb->i.x - s1->i.x, + x21 = s2->i.x - s1->i.x, + yb1 = sb->i.y - s1->i.y, + y21 = s2->i.y - s1->i.y, + xba = sb->i.x - sa->i.x, + x2a = s2->i.x - sa->i.x, + yba = sb->i.y - sa->i.y, + y2a = s2->i.y - sa->i.y; + register double + cosb12 = double(xb1*x21 + yb1*y21), + cosba2 = double(xba*x2a + yba*y2a) , + sinb12 = double(det2), + sinba2 = double(t2->det); + + + // angle b12 > angle ba2 => cotg(angle b12) < cotg(angle ba2) + OnSwap = ((double) cosb12 * (double) sinba2) < ((double) cosba2 * (double) sinb12); + break; + } + else + { + // critere de Delaunay anisotrope + Real8 som; + I2 AB=(I2) *sb - (I2) *sa; + I2 MAB2=((I2) *sb + (I2) *sa); + R2 MAB(MAB2.x*0.5,MAB2.y*0.5); + I2 A1=(I2) *s1 - (I2) *sa; + I2 D = (I2) * s1 - (I2) * sb ; + R2 S2(s2->i.x,s2->i.y); + R2 S1(s1->i.x,s1->i.y); + { + Metric M=s1->m; + R2 ABo = M.Orthogonal(AB); + R2 A1o = M.Orthogonal(A1); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + if (Abs(d) > dd*1.e-3) { + R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + som = M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + + } + { + Metric M=s2->m; + R2 ABo = M.Orthogonal(AB); + R2 A1o = M.Orthogonal(A1); + // (A+B)+ x ABo = (S1+B)/2+ y A1 + // ABo x - A1o y = (S1+B)/2-(A+B)/2 = (S1-B)/2 = D/2 + double dd = Abs(ABo.x*A1o.y)+Abs(ABo.y*A1o.x); + double d = (ABo.x*A1o.y - ABo.y*A1o.x)*2; // because D/2 + if(Abs(d) > dd*1.e-3) { + R2 C(MAB+ABo*((D.x*A1o.y - D.y*A1o.x)/d)); + som += M(C - S2)/M(C - S1); + } else + {kopt=1;continue;} + } + OnSwap = som < 2; + break; + } + + } // OnSwap + } // (! OnSwap &&(det1 > 0) && (det2 > 0) ) + } + if( OnSwap ) + bamg::swap(t1,a1,t2,a2,s1,s2,det1,det2); + else { + NbUnSwap ++; + t1->SetMarkUnSwap(a1); + } + return OnSwap; } /*}}}1*/ - /*FUNCTION Triangles::ShowHistogram{{{1*/ - void Triangles::ShowHistogram() const { - - const Int4 kmax=10; - const Real8 llmin = 0.5,llmax=2; - const Real8 lmin=log(llmin),lmax=log(llmax),delta= kmax/(lmax-lmin); - Int4 histo[kmax+1]; - Int4 i,it,k, nbedges =0; - for (i=0;i<=kmax;i++) histo[i]=0; - for (it=0;itlink || Number(ta) >= it) - { - Vertex & vP = triangles[it][VerticesOfTriangularEdge[j][0]]; - Vertex & vQ = triangles[it][VerticesOfTriangularEdge[j][1]]; - if ( !& vP || !&vQ) continue; - R2 PQ = vQ.r - vP.r; - Real8 l = log(LengthInterpole(vP,vQ,PQ)); - nbedges++; - k = (int) ((l - lmin)*delta); - k = Min(Max(k,0L),kmax); - histo[k]++; - } - } - } - printf(" --- Histogram of the unit mesh, nb of edges = %i\n",nbedges); - printf(" length of edge in | %% of edge | Nb of edges \n"); - printf(" --------------------+-------------+-------------\n"); - for (i=0;i<=kmax;i++){ - if (i==0) printf(" %10i",0); - else printf(" %10g",exp(lmin+i/delta)); - if (i==kmax) printf(" +inf "); - else printf(" %10g",exp(lmin+(i+1)/delta)); - printf("| %10g |\n",((long) ((10000.0 * histo[i])/ nbedges))/100.0); - printf(" %g\n",histo[i]); - } - printf(" --------------------+-------------+-------------\n"); +/*FUNCTION Triangles::UnCrack{{{1*/ +int Triangles::UnCrack() { + if (NbCrackedEdges!=0 && NbCrackedVertices<=0);{ + throw ErrorException(__FUNCT__,exprintf("NbCrackedEdges ==0 || NbCrackedVertices >0")); } - /*}}}1*/ - /*FUNCTION Triangles::Crack{{{1*/ - int Triangles::Crack() { - if (NbCrackedEdges!=0 && NbCrackedVertices<=0);{ - throw ErrorException(__FUNCT__,exprintf("NbCrackedEdges!=0 && NbCrackedVertices<=0")); - } - for (int i=0;i0")); - } - for (int i=0;iCracked()) k++; - if( k==0) return 0; - CurrentTh = this; - printf(" number of Cracked Edges = %i\n",k); - NbCrackedEdges =k; - CrackedEdges = new CrackedEdge[k]; - // new edge - Edge * e = new Edge[ nbe + k]; - - // copy - for (i=0;iCracked()) - { - e[nbe] = e[i]; - // return the edge - e[nbe].v[0] = e[i].v[1]; - e[nbe].v[1] = e[i].v[0]; - e[nbe].on = e[i].on->link ; // fqux - CrackedEdges[k++]=CrackedEdge(edges,i,nbe); - nbe++; - } - ReMakeTriangleContainingTheVertex() ; - // - int nbcrakev =0; - Vertex *vlast = vertices + nbv; - Vertex *vend = vertices + nbvx; // end of array - for (int iv=0;iv=3)){ - throw ErrorException(__FUNCT__,exprintf("!tbegin || (i<0) || (i>=3)")); - } - // turn around the vertex v - TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge - int k=0; - do { - int kv = VerticesOfTriangularEdge[ta][1]; - k++; - Triangle * tt (ta); - if ( ta.Cracked() ) - { - TriangleAdjacent tta=(ta.Adj()); - if (!tta.Cracked()){ - throw ErrorException(__FUNCT__,exprintf("!tta.Cracked()")); - } - if ( kk == 0) tbegin=ta,kkk=0; // begin by a cracked edge => restart - if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use - kk++;// number of cracked edge view - } - if ( tt->link ) { // if good triangles store the value - int it = Number(tt); - if (it>=nt){ - throw ErrorException(__FUNCT__,exprintf("(it>=nt)")); - } - (*tt)(kv)= vv; // Change the vertex of triangle - if(vvReferenceNumber=iv;} // copy the vertex value + store the old vertex number in ref - // tt->SetTriangleContainingTheVertex(); - kkk++; - } else if (kk) { // crack + boundary - if ( kkk ) { kc =1;vv = vlast++; kkk = 0; } // new vertex if use - } - - ta = Next(ta).Adj(); - } while ( (tbegin != ta)); - if (!k){ - throw ErrorException(__FUNCT__,exprintf("!k")); - } - if (kc) nbcrakev++; - } - - if ( nbcrakev ) - for (int iec =0;iec < NbCrackedEdges; iec ++) - CrackedEdges[iec].Set(); - - // set the ref - NbCrackedVertices = nbcrakev; - // int nbvo = nbv; - nbv = vlast - vertices; - int nbnewv = nbv - nbv; // nb of new vrtices - if (nbcrakev && verbosity > 1 ) printf(" number of Cracked vertices = %i, number of created vertices = %i\n",nbcrakev,nbnewv); - // all the new vertices are on geometry - if (nbnewv) - { // - Int4 n = nbnewv+NbVerticesOnGeomVertex; - Int4 i,j,k; - VertexOnGeom * vog = new VertexOnGeom[n]; - for ( i =0; i= LastOld) - { // a new vertex - Int4 old = v->ReferenceNumber ; // the old same vertex - Int4 i = ( v - LastOld); - // if the old is on vertex => warning - // else the old is on edge => ok - vog[i] = vog[old]; - // vog[i].mv = v; - //g[i].ge = ; - //og[i].abcisse = ; - } - - } - } - - NbVerticesOnGeomVertex = n; - } - SetVertexFieldOn(); - - if (vlast >= vend) { - throw ErrorException(__FUNCT__,exprintf("Not enougth vertices: to crack the mesh we need %i vertices",nbv)); - } - CurrentTh = CurrentThOld; - return NbCrackedVertices; - } - /*}}}1*/ - /*FUNCTION Triangles::FindTriangleContening{{{1*/ - Triangle * Triangles::FindTriangleContening(const I2 & B,Icoor2 dete[3], Triangle *tstart) const { - Triangle * t=0; - int j,jp,jn,jj; - if (tstart) t=tstart; - else { - if (!quadtree){ - throw ErrorException(__FUNCT__,exprintf("!quadtree")); - } - Vertex *a = quadtree->NearestVertex(B.x,B.y) ; - - if (! a || !a->t ) { - if (a) { - printf("TriangleConteningTheVertex vertex number %i, another call to ReMakeTriangleContainingTheVertex was required\n", Number(a)); - } - throw ErrorException(__FUNCT__,exprintf("problem in Triangles::FindTriangleContening")); - } - if (a=vertices+nbv){ - throw ErrorException(__FUNCT__,exprintf("a=vertices+nbv")); - } - t = a->t; - if (t=triangles+nbt){ - throw ErrorException(__FUNCT__,exprintf("t=triangles+nbt")); - } - } - Icoor2 detop ; - int kkkk =0; // number of test triangle - - while ( t->det < 0){ // the initial triangles is outside - int k0=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; - if (k0<0){ // k0 the NULL vertex - throw ErrorException(__FUNCT__,exprintf("k0<0")); - } - 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 - return t; - t = t->TriangleAdj(OppositeEdge[k0]); - kkkk++; - if (kkkk>=2){ - throw ErrorException(__FUNCT__,exprintf("kkkk>=2")); - } - } - - jj=0; - detop = det(*(*t)(VerticesOfTriangularEdge[jj][0]),*(*t)(VerticesOfTriangularEdge[jj][1]),B); - - while(t->det > 0 ) { - kkkk++; - if (kkkk>=2000){ - throw ErrorException(__FUNCT__,exprintf("kkkk>=2000")); - } - j= OppositeVertex[jj]; - dete[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 - 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; - // 0 => ok - // 1 => go in way 1 - // 2 => two way go in way 1 or 2 randomly - - if (k==0) break; - if (k == 2 && BinaryRand()) Exchange(ii[0],ii[1]); - if ( k>=3){ - throw ErrorException(__FUNCT__,exprintf("k>=3")); - } - TriangleAdjacent t1 = t->Adj(jj=ii[0]); - if ((t1.det() < 0 ) && (k == 2)) - t1 = t->Adj(jj=ii[1]); - t=t1; - j=t1;// for optimisation we now the -det[OppositeVertex[j]]; - detop = -dete[OppositeVertex[jj]]; - jj = j; - } - - if (t->det<0) // outside triangle - dete[0]=dete[1]=dete[2]=-1,dete[OppositeVertex[jj]]=detop; - // NbOfTriangleSearchFind += kkkk; - return t; - } - /*}}}1*/ - /*FUNCTION Triangles::IntersectGeomMetric{{{1*/ - void Triangles::IntersectGeomMetric(const Real8 err=1,const int iso=0){ - long int verbosity=0; - Real8 ss[2]={0.00001,0.99999}; - Real8 errC = 2*sqrt(2*err); - Real8 hmax = Gh.MaximalHmax(); - Real8 hmin = Gh.MinimalHmin(); - Real8 maxaniso = 1e6; - if (hmax<=0){ - throw ErrorException(__FUNCT__,exprintf("hmax<=0")); - } - SetVertexFieldOn(); - if (errC > 1) errC = 1; - for (Int4 i=0;iR1tg(s,tg); - Real8 ht = hmax; - if (R1>1.0e-20) { // err relative to the length of the edge - ht = Min(Max(errC/R1,hmin),hmax); - } - Real8 hn = iso? ht : Min(hmax,ht*maxaniso); - if (ht<=0 || hn<=0){ - throw ErrorException(__FUNCT__,exprintf("ht<=0 || hn<=0")); - } - MatVVP2x2 Vp(1/(ht*ht),1/(hn*hn),tg); - Metric MVp(Vp); - edges[i][j].m.IntersectWith(MVp); - } - - } - // the problem is for the vertex on vertex - } - /*}}}1*/ - /*FUNCTION Triangles::BoundAnisotropy{{{1*/ - void Triangles::BoundAnisotropy(Real8 anisomax,Real8 hminaniso) { - long int verbosity=0; - - double lminaniso = 1/ (Max(hminaniso*hminaniso,1e-100)); - if (verbosity > 1) printf(" BoundAnisotropy by %g\n",anisomax); - Real8 h1=1.e30,h2=1e-30,rx=0; - Real8 coef = 1./(anisomax*anisomax); - Real8 hn1=1.e30,hn2=1e-30,rnx =1.e-30; - for (Int4 i=0;i2){ - printf(" input: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5),pow(h1,-0.5),pow(rx,0.5)); - printf(" output: Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); - } - } - /*}}}1*/ - /*FUNCTION Triangles::IntersectConsMetric{{{1*/ - void Triangles::IntersectConsMetric(const double * s,const Int4 nbsol,const int * typsols, - const Real8 hmin1,const Real8 hmax1,const Real8 coef, - const Real8 anisomax ,const Real8 CutOff,const int NbJacobi, - const int DoNormalisation,const double power,const int choice) - { // the array of solution s is store - // sol0,sol1,...,soln on vertex 0 - // sol0,sol1,...,soln on vertex 1 - // etc. - // choise = 0 => H is computed with green formule - // otherwise => H is computed from P2 on 4T - const int dim = 2; - - long int verbosity=0; - - int sizeoftype[] = { 1, dim ,dim * (dim+1) / 2, dim * dim } ; - - // computation of the nb of field - Int4 ntmp = 0; - if (typsols) - { - for (Int4 i=0;i1e-30; - Real8 hmin = Max(hmin1,MinimalHmin()); - Real8 hmax = Min(hmax1,MaximalHmax()); - Real8 coef2 = 1/(coef*coef); - - if(verbosity>1) { - printf(" Construction of Metric: number of field: %i (nbt=%i, nbv=%i)\n",n,nbt,nbv); - printf(" coef = %g\n",coef); - printf(" hmin = %g hmax = %g\n",hmin,hmax); - printf(" anisomax = %g nb Jacobi = %i, power = %i\n",anisomax,NbJacobi,power); - if (RelativeMetric) printf(" RelativeErr with CutOff= %g\n",CutOff); - else printf(" Absolute error\n"); - } - double *ss=(double*)s; - - double sA,sB,sC; - - Real8 *detT = new Real8[nbt]; - Real8 *Mmass= new Real8[nbv]; - Real8 *Mmassxx= new Real8[nbv]; - Real8 *dxdx= new Real8[nbv]; - Real8 *dxdy= new Real8[nbv]; - Real8 *dydy= new Real8[nbv]; - Real8 *workT= new Real8[nbt]; - Real8 *workV= new Real8[nbv]; - int *OnBoundary = new int[nbv]; - for (iv=0;ivlink) // no adj triangle => edge on boundary - OnBoundary[Number(t[VerticesOfTriangularEdge[j][0]])]=1, - OnBoundary[Number(t[VerticesOfTriangularEdge[j][1]])]=1, - nbb++; - } - - workT[i] = nbb; - Mmass[iA] += dett; - Mmass[iB] += dett; - Mmass[iC] += dett; - - if((nbb==0)|| !choice) - { - Mmassxx[iA] += dett; - Mmassxx[iB] += dett; - Mmassxx[iC] += dett; - } - } - else - workT[i]=-1; - - for (Int4 nusol=0;nusol2) printf(" Solution %i, Min = %g, Max = %g, Delta = %g, cnorm = %g, number of fields = %i\n",nusol,smin,smax,sdelta,cnorm,nbfield); - - if ( sdelta < 1.0e-10*Max(absmax,1e-20) && (nbfield ==1)) { - if (verbosity>2) printf(" Solution %i is constant, skipping...\n"); - continue; - } - - double *sf = ss; - for (Int4 nufield=0;nufieldlink) - { - Vertex &v = *ta.OppositeVertex(); - R2 V = v; - Int4 iV = Number(v); - Real8 lA = bamg::Area2(V,B,C)/dd; - Real8 lB = bamg::Area2(A,V,C)/dd; - Real8 lC = bamg::Area2(A,B,V)/dd; - taa[0][j] = lB*lC; - taa[1][j] = lC*lA; - taa[2][j] = lA*lB; - lla = lA,llb=lB,llc=lC,llf=ss[iV*n] ; - - bb[j] = ss[iV*n] - ( sA*lA + sB*lB + sC*lC ) ; - } - else - { - nbb++; - taa[0][j]=0; - taa[1][j]=0; - taa[2][j]=0; - taa[j][j]=1; - bb[j]=0; - } - } - - // resolution of 3x3 lineaire system transpose - Real8 det33 = det3x3(taa[0],taa[1],taa[2]); - Real8 cBC = det3x3(bb,taa[1],taa[2]); - Real8 cCA = det3x3(taa[0],bb,taa[2]); - Real8 cAB = det3x3(taa[0],taa[1],bb); - - if (!det33){ - throw ErrorException(__FUNCT__,exprintf("!det33")); - } - // computation of the gradient in the element - - // H( li*lj) = grad li grad lj + grad lj grad lj - // grad li = njk / detT ; with i j k =(A,B,C) - Real8 Hxx = cAB * ( nBC.x*nCA.x) + cBC * ( nCA.x*nAB.x) + cCA * (nAB.x*nBC.x); - Real8 Hyy = cAB * ( nBC.y*nCA.y) + cBC * ( nCA.y*nAB.y) + cCA * (nAB.y*nBC.y); - Real8 Hxy = cAB * ( nBC.y*nCA.x) + cBC * ( nCA.y*nAB.x) + cCA * (nAB.y*nBC.x) - + cAB * ( nBC.x*nCA.y) + cBC * ( nCA.x*nAB.y) + cCA * (nAB.x*nBC.y); - Real8 coef = 1.0/(3*dd*det33); - Real8 coef2 = 2*coef; - Hxx *= coef2; - Hyy *= coef2; - Hxy *= coef2; - if(nbb==0) - { - dxdx[iA] += Hxx; - dydy[iA] += Hyy; - dxdy[iA] += Hxy; - - dxdx[iB] += Hxx; - dydy[iB] += Hyy; - dxdy[iB] += Hxy; - - dxdx[iC] += Hxx; - dydy[iC] += Hyy; - dxdy[iC] += Hxy; - } - - } - else - { - - // if edge on boundary no contribution => normal = 0 - if ( ! tBC || ! tBC->link ) nBC = O; - if ( ! tCA || ! tCA->link ) nCA = O; - if ( ! tAB || ! tAB->link ) nAB = O; - - // remark we forgot a 1/2 because - // $\\int_{edge} w_i = 1/2 $ if $i$ is in edge - // 0 if not - // if we don't take the boundary - // dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; - - dxdx[iA] += ( nCA.x + nAB.x ) *Grads.x; - dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x; - dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x; - - // warning optimization (1) the divide by 2 is done on the metrix construction - dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ; - dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ; - dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ; - - dydy[iA] += ( nCA.y + nAB.y ) *Grads.y; - dydy[iB] += ( nAB.y + nBC.y ) *Grads.y; - dydy[iC] += ( nBC.y + nCA.y ) *Grads.y; - } - - } // for real all triangles - Int4 kk=0; - for ( iv=0,k=0 ; iv0) - { - dxdx[iv] /= 2*Mmassxx[iv]; - // warning optimization (1) on term dxdy[iv]*ci/2 - dxdy[iv] /= 4*Mmassxx[iv]; - dydy[iv] /= 2*Mmassxx[iv]; - // Compute the matrix with abs(eigen value) - Metric M(dxdx[iv], dxdy[iv], dydy[iv]); - MatVVP2x2 Vp(M); - Vp.Abs(); - M = Vp; - dxdx[iv] = M.a11; - dxdy[iv] = M.a21; - dydy[iv] = M.a22; - } - else kk++; - - - // correction of second derivate - // by a laplacien - - Real8 *d2[3] = { dxdx, dxdy, dydy}; - Real8 *dd; - for (int xy = 0;xy<3;xy++) - { - dd = d2[xy]; - // do leat 2 iteration for boundary problem - for (int ijacobi=0;ijacobi0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.); - workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc; - } - for (iv=0;iv2) { - printf(" Field %i of solution %i\n",nufield,nusol); - printf(" before bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(h2,-0.5), pow(h1,-0.5), pow(rx,0.5)); - printf(" after bounding : Hmin = %g, Hmax = %g, factor of anisotropy max = %g\n",pow(hn2,-0.5),pow(hn1,-0.5),pow(rnx,0.5)); - } - } // end of for all field - }// end for all solution - - delete [] detT; - delete [] Mmass; - delete [] dxdx; - delete [] dxdy; - delete [] dydy; - delete [] workT; - delete [] workV; - delete [] Mmassxx; - delete [] OnBoundary; - - } - /*}}}1*/ - /*FUNCTION Triangles::MaxSubDivision{{{1*/ - void Triangles::MaxSubDivision(Real8 maxsubdiv) { - long int verbosity=0; - - const Real8 maxsubdiv2 = maxsubdiv*maxsubdiv; - if(verbosity>1) printf(" Limit the subdivision of a edges in the new mesh by %g\n",maxsubdiv); - // for all the edges - // if the len of the edge is to long - Int4 it,nbchange=0; - Real8 lmax=0; - for (it=0;it maxsubdiv2) - { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M - Real8 lc = M(AC,AC); - D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; - D2xD2 Rt1(Rt.inv()); - D2xD2 D(maxsubdiv2,0,0,lc); - D2xD2 MM = Rt1*D*Rt1.t(); - v0.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y); - nbchange++; - } - M = v1; - l = M(AB,AB); - lmax = Max(lmax,l); - if(l> maxsubdiv2) - { R2 AC = M.Orthogonal(AB);// the ortogonal vector of AB in M - Real8 lc = M(AC,AC); - D2xD2 Rt(AB,AC);// Rt.x = AB , Rt.y = AC; - D2xD2 Rt1(Rt.inv()); - D2xD2 D(maxsubdiv2,0,0,lc); - D2xD2 MM = Rt1*D*Rt1.t(); - v1.m = M = MetricAnIso(MM.x.x,MM.y.x,MM.y.y); - nbchange++; - } - - - } - } - } - if(verbosity>3){ - printf(" number of metric changes = %i, maximum number of subdivision of a edges before change = %g\n",nbchange,pow(lmax,0.5)); - } - } - /*}}}1*/ - /*FUNCTION Triangles::SmoothMetric{{{1*/ - void Triangles::SmoothMetric(Real8 raisonmax) { - long int verbosity=0; - - if(raisonmax<1.1) return; - if(verbosity > 1) printf(" Triangles::SmoothMetric raisonmax = %g\n",raisonmax); - ReMakeTriangleContainingTheVertex(); - Int4 i,j,kch,kk,ip; - Int4 *first_np_or_next_t0 = new Int4[nbv]; - Int4 *first_np_or_next_t1 = new Int4[nbv]; - Int4 Head0 =0,Head1=-1; - Real8 logseuil= log(raisonmax); - - for(i=0;i=0&& kk++<100) { - kch=0; - for (i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) { - // pour tous les triangles autour du sommet s - register Triangle* t= vertices[i].t; - if (!t){ - throw ErrorException(__FUNCT__,exprintf("!t")); - } - Vertex & vi = vertices[i]; - TriangleAdjacent ta(t,EdgesVertexTriangle[vertices[i].vint][0]); - Vertex *pvj0 = ta.EdgeVertex(0); - while (1) { - ta=Previous(Adj(ta)); - if (vertices+i != ta.EdgeVertex(1)){ - throw ErrorException(__FUNCT__,exprintf("vertices+i != ta.EdgeVertex(1)")); - } - Vertex & vj = *(ta.EdgeVertex(0)); - if ( &vj ) { - j= &vj-vertices; - if (j<0 || j >= nbv){ - throw ErrorException(__FUNCT__,exprintf("j<0 || j >= nbv")); - } - R2 Aij = (R2) vj - (R2) vi; - Real8 ll = Norme2(Aij); - if (0) { - Real8 hi = ll/vi.m(Aij); - Real8 hj = ll/vj.m(Aij); - if(hi < hj) - { - Real8 dh=(hj-hi)/ll; - if (dh>logseuil) { - vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi)); - if(first_np_or_next_t1[j]<0) - kch++,first_np_or_next_t1[j]=Head1,Head1=j; - } - } - } - else - { - Real8 li = vi.m(Aij); - if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) ) - if(first_np_or_next_t1[j]<0) // if the metrix change - kch++,first_np_or_next_t1[j]=Head1,Head1=j; - } - } - if ( &vj == pvj0 ) break; - } - } - Head0 = Head1; - Head1 = -1; - Exchange(first_np_or_next_t0,first_np_or_next_t1); - } - if(verbosity>2) printf(" number of iterations = %i\n",kch); - delete [] first_np_or_next_t0; - delete [] first_np_or_next_t1; - } - /*}}}1*/ - /*FUNCTION Triangles::NearestVertex{{{1*/ - Vertex * Triangles::NearestVertex(Icoor1 i,Icoor1 j) { - return quadtree->NearestVertex(i,j); - } - /*}}}1*/ - /*FUNCTION Triangles::PreInit{{{1*/ - void Triangles::PreInit(Int4 inbvx,char *fname) { - long int verbosity=0; - - srand(19999999); - OnDisk =0; - NbRef=0; - // allocGeometry=0; - identity=0; - NbOfTriangleSearchFind =0; - NbOfSwapTriangle =0; - nbiv=0; - nbv=0; - nbvx=inbvx; - nbt=0; - NbOfQuad = 0; - nbtx=2*inbvx-2; - NbSubDomains=0; - NbVertexOnBThVertex=0; - NbVertexOnBThEdge=0; - VertexOnBThVertex=0; - VertexOnBThEdge=0; - - NbCrackedVertices=0; - NbCrackedEdges =0; - CrackedEdges =0; - nbe = 0; - name = fname ; - - if (inbvx) { - vertices=new Vertex[nbvx]; - if (!vertices){ - throw ErrorException(__FUNCT__,exprintf("!vertices")); - } - ordre=new (Vertex* [nbvx]); - if (!ordre){ - throw ErrorException(__FUNCT__,exprintf("!ordre")); - } - triangles=new Triangle[nbtx]; - if (!triangles){ - throw ErrorException(__FUNCT__,exprintf("!triangles")); - } - } - else { - vertices=0; - ordre=0; - triangles=0; - nbtx=0; - } - if ( name || inbvx) { - time_t timer =time(0); - char buf[70]; - strftime(buf ,70,", Date: %y/%m/%d %H:%M %Ss",localtime(&timer)); - counter++; - char countbuf[30]; - sprintf(countbuf,"%d",counter); - int lg =0 ; - if (&BTh != this && BTh.name) - lg = strlen(BTh.name)+4; - identity = new char[ lg + strlen(buf) + strlen(countbuf)+ 2 + 10 + ( Gh.name ? strlen(Gh.name) + 4 : 0)]; - identity[0]=0; - if (lg) - strcat(strcat(strcat(identity,"B="),BTh.name),", "); - - if (Gh.name) - strcat(strcat(identity,"G="),Gh.name); - strcat(strcat(identity,";"),countbuf); - strcat(identity,buf); - } - - quadtree=0; - edges=0; - VerticesOnGeomVertex=0; - VerticesOnGeomEdge=0; - NbVerticesOnGeomVertex=0; - NbVerticesOnGeomEdge=0; - subdomains=0; - NbSubDomains=0; - } - /*}}}1*/ + for (int i=0;i