Index: /issm/trunk/src/c/Bamgx/Mesh2.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2842) +++ /issm/trunk/src/c/Bamgx/Mesh2.cpp (revision 2843) @@ -1,31 +1,2 @@ -// -*- Mode : c++ -*- -// -// SUMMARY : -// USAGE : -// ORG : -// AUTHOR : Frederic Hecht -// E-MAIL : hecht@ann.jussieu.fr -// - -/* - - This file is part of Freefem++ - - Freefem++ is free software; you can redistribute it and/or modify - it under the terms of the GNU Lesser General Public License as published by - the Free Software Foundation; either version 2.1 of the License, or - (at your option) any later version. - - Freefem++ is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public License - along with Freefem++; if not, write to the Free Software - Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - */ -extern bool withrgraphique; - #include "../shared/shared.h" #include @@ -35,36 +6,18 @@ #include -using namespace std; - #include "Mesh2.h" #include "QuadTree.h" #include "SetOfE4.h" +using namespace std; namespace bamg { int Triangles::counter = 0; - - Triangles * CurrentTh =0; - + Triangles* CurrentTh=0; int hinterpole=1; - int ForDebugging = 0; const Direction NoDirOfSearch = Direction(); -#ifndef NDEBUG - inline void MyAssert(int i,char*ex,char * file,long line) - { - if( i) { - cerr << "Error Assert:" << ex << " in " << file << " line: " << line << endl; -#ifdef NOTFREEFEM - exit(1); -#else - throw(ErrorExec("exit",1000)); -#endif - } - } -#endif - - Int4 AGoodNumberPrimeWith(Int4 n) - { + + Int4 AGoodNumberPrimeWith(Int4 n){ const Int4 BigPrimeNumber[] ={ 567890359L, 567890431L, 567890437L, 567890461L, 567890471L, @@ -83,16 +36,10 @@ } - class Triangles; void MeshError(int Err,Triangles *Th){ cerr << " Fatal error in the meshgenerator " << Err << endl ; -#ifdef NOTFREEFEM exit(1); -#else - throw(ErrorMesh("Bamg",Err,Th)); -#endif } - ostream& operator <<(ostream& f, const Triangle & ta) - { + ostream& operator <<(ostream& f, const Triangle & ta) { if(CurrentTh) f << "[" << CurrentTh->Number(ta) << "::" @@ -115,653 +62,4 @@ return f;} - void swap(Triangle *t1,Int1 a1, - Triangle *t2,Int1 a2, - Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2) - { // swap - // -------------------------------------------------------------- - // Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles - // - // sb sb - // / | \ / \ ! - // as1/ | \ /a2 \ ! - // / | \ / t2 \ ! - // s1 /t1 | t2 \s2 --> s1 /___as2___\s2 ! - // \ a1|a2 / \ as1 / - // \ | / \ t1 / - // \ | / as2 \ a1/ - // \ | / \ / - // sa sa - // ------------------------------------------------------------- - int as1 = NextEdge[a1]; - int as2 = NextEdge[a2]; - int ap1 = PreviousEdge[a1]; - int ap2 = PreviousEdge[a2]; - (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb - (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa - // mise a jour des 2 adjacences externes - TriangleAdjacent taas1 = t1->Adj(as1), - taas2 = t2->Adj(as2), - tas1(t1,as1), tas2(t2,as2), - ta1(t1,a1),ta2(t2,a2); - // externe haut gauche - taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap()); - // externe bas droite - taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap()); - // remove the Mark UnMarkSwap - t1->SetUnMarkUnSwap(ap1); - t2->SetUnMarkUnSwap(ap2); - // interne - tas1.SetAdj2(tas2); - - t1->det = det1; - t2->det = det2; - - t1->SetTriangleContainingTheVertex(); - t2->SetTriangleContainingTheVertex(); - } // end swap - - - - - - Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside) - { - CurrentTh=&Th; - assert(&Th); - I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); - Icoor2 dete[3]; - Triangle & tb = *Th.FindTriangleContening(I,dete); - - if (tb.link) - { // internal point in a true triangles - a[0]= (Real8) dete[0]/ tb.det; - a[1]= (Real8) dete[1] / tb.det; - a[2] = (Real8) dete[2] / tb.det; - inside = 1; - return Th.Number(tb); - } - else - { - inside = 0; - double aa,bb; - TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb); - int k = ta; - Triangle * tc = ta; - if (!tc->link) - { ta = ta.Adj(); - tc=ta; - k = ta; - Exchange(aa,bb); - assert(tc->link); - } - a[VerticesOfTriangularEdge[k][0]] = aa; - a[VerticesOfTriangularEdge[k][1]] = bb; - a[OppositeVertex[k]] = 1- aa -bb; - return Th.Number(tc); - } - } - - - TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) { - // - // cout << " - "; - int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; - int dir=0; - assert(k>=0); - int kkk=0; - Icoor2 IJ_IA,IJ_AJ; - TriangleAdjacent edge(t,OppositeEdge[k]); - for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) - { - - assert(kkk++<1000); - Vertex &vI = *edge.EdgeVertex(0); - Vertex &vJ = *edge.EdgeVertex(1); - I2 I=vI, J=vJ, IJ= J-I; - IJ_IA = (IJ ,(A-I)); - // cout << A << vI.i << vJ.i << edge << " " << IJ_IA << " dir " << dir <0) {a=1;b=0;return edge;}// change of signe => I - else {dir=-1; - continue;}};// go in direction i - IJ_AJ = (IJ ,(J-A)); - if (IJ_AJ<0) { - if(dir<0) {a=0;b=1;return edge;} - else {dir = 1; - continue;}}// go in direction j - double IJ2 = IJ_IA + IJ_AJ; - assert(IJ2); - a= IJ_AJ/IJ2; - b= IJ_IA/IJ2; - // cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl; - return edge; - } - } - - TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) - { - // walk around the vertex - // version 2 for remove the probleme if we fill the hole - //int bug=1; - // Triangle *torigine = t; - // restart: - // int dir=0; - assert(t->link == 0); - // to have a starting edges - // try the 3 edge bourna-- in case of internal hole - // and choice the best - // - // - // the probleme is in case of the fine and long internal hole - // for exemple neart the training edge of a wing - // - Vertex * s=0,*s1=0, *s0=0; - Icoor2 imax = MaxICoor22; - Icoor2 l0 = imax,l1 = imax; - double dd2 = imax;// infinity - TriangleAdjacent er; - int cas=-2; - for (int j=0;j<3;j++) - { - TriangleAdjacent ta=t->FindBoundaryEdge(j); - if (! (Triangle *) ta) continue; - s0 = ta.EdgeVertex(0); - s1 = ta.EdgeVertex(1); - I2 A = * s0; - I2 B = *ta.EdgeVertex(1); - I2 AB = B-A,AC=C-A,BC=B-C; - Icoor2 ACAC = (AC,AC), BCBC = (BC,BC); - Icoor2 AB2 = Norme2_2(AB); // ||AB||^2 - Icoor2 ABAC = (AB,AC); // AB.AC| - - double d2; - if ( ABAC < 0 ) // DIST A - { - if ( (d2=(double) ACAC) < dd2) - { - // cout << " A " << d2 << " " << dd2; - er = ta; - l0 = ACAC; - l1 = BCBC; - cas = 0; - s = s0; - } - } - else if (ABAC > AB2) // DIST B - { - if ( (d2=(double) BCBC) < dd2) - { - // cout << " B " << d2 << " " << dd2; - dd2 = d2; - er = Adj(ta); // other direction - l0 = BCBC; - l1 = ACAC; - cas = 1; - s = s1; - } - } - else // DIST AB - { - - double det_2 = (double) Det(AB,AC); - det_2 *= det_2; // square of area*2 of triangle ABC - d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC - // cout << " AB " << d2 << " " << dd2 - // << " " << CurrentTh->Number(ta.EdgeVertex(0)) - // << " " << CurrentTh->Number(ta.EdgeVertex(1)) << " " ; - - if (d2 < dd2) - { - dd2 = d2; - er = ta; - l0 = (AC,AC); - l1 = (BC,BC); - s = 0; - cas = -1; - // cout << " ABAC " << ABAC << " ABAC " << ABAC - // << " AB2 " << AB2 << endl; - b = ((double) ABAC/(double) AB2); - a = 1 - b; - } - } - } - assert(cas !=-2); - // l1 = ||C s1|| , l0 = ||C s0|| - // where s0,s1 are the vertex of the edge er - - if ( s) - { - t=er; - TriangleAdjacent edge(er); - - int kkk=0; - int linkp = t->link == 0; - - Triangle * tt=t=edge=Adj(Previous(edge)); - // cout << CurrentTh->Number(t) << " " << linkp << endl; - do { // loop around vertex s - - assert(edge.EdgeVertex(0)==s && kkk++<10000); - - int link = tt->link == 0; - // cout << CurrentTh->Number(tt) << " " << link << " " << CurrentTh->Number(s) - // << " " << CurrentTh->Number(er.EdgeVertex(0)) - // << " " << CurrentTh->Number(er.EdgeVertex(1)) - // << " " << CurrentTh->Number(edge.EdgeVertex(0)) - // << " " << CurrentTh->Number(edge.EdgeVertex(1)) - // << endl; - if ((link + linkp) == 1) - { // a boundary edge - Vertex * st = edge.EdgeVertex(1); - I2 I=*st; - Icoor2 ll = Norme2_2 (C-I); - if (ll < l1) { // the other vertex is neart - s1=st; - l1=ll; - er = edge; - if(llNumber(er.EdgeVertex(0)) - // << " " << CurrentTh->Number(er.EdgeVertex(1)) ; - if (aa<0) a=1,b=0; - else if(bb<0) a=0,b=1; - else - { - a = bb/(aa+bb); - b = aa/(aa+bb); - } - } - - // cout <<" return= " << CurrentTh->Number(er.EdgeVertex(0)) << " " - // << CurrentTh->Number(er.EdgeVertex(1)) << " " << a - // << " " << b <<" " << l0 << " " <det>=0) - ifirst = ilast;} - else {// not optimisation - init(); - t=tbegin = Bh.FindTriangleContening(a,deta); - if( t->det>=0) - ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det); - else - {// find the nearest boundary edge of the vertex A - // find a edge or such normal projection a the edge IJ is on the edge - // <=> IJ.IA >=0 && IJ.AJ >=0 - ilast=ifirst; - double ba,bb; - TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); - Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); - NewItem(A,Metric(ba,v0,bb,v1)); - t=edge; - // test if the point b is in the same side - if (det(v0.i,v1.i,b)>=0) { - //cout << " All the edge " << A << B << endl; - TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); - Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); - NewItem(A,Metric(ba,v0,bb,v1)); - return; - } - } // find the nearest boundary edge of the vertex A - } // end not optimisation - if (t->det<0) { // outside departure - while (t->det <0) { // intersection boundary edge and a,b, - k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; - assert(k>=0); - ocut = OppositeEdge[k]; - i=VerticesOfTriangularEdge[ocut][0]; - j=VerticesOfTriangularEdge[ocut][1]; - vi=(*t)[i]; - vj=(*t)[j]; - deti = bamg::det(a,b,vi); - detj = bamg::det(a,b,vj); - if (deti>0) // go to i direction on gamma - ocut = PreviousEdge[ocut]; - else if (detj<=0) // go to j direction on gamma - ocut = NextEdge[ocut]; - TriangleAdjacent tadj =t->Adj(ocut); - t = tadj; - iedge= tadj; - if (t == tbegin) { // - double ba,bb; - long int verbosity=2; - if (verbosity>7) - cout << " SplitEdge: All the edge " << A << B << nbegin << det(vi,vj,b) - << " deti= " << deti << " detj=" <det <0) - // theoriticaly we have: deti =<0 and detj>0 - - // computation of barycentric coor - // test if the point b is on size on t - // we revert vi,vj because vi,vj is def in Adj triangle - if ( det(vi,vj,b)>=0) { - if (verbosity>7) - cout << " SplitEdge: all AB outside " << A << B << endl; - t=tbegin; - Real8 ba,bb; - TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb); - NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1))); - return; - } - else - { - k = OppositeVertex[iedge]; - i=VerticesOfTriangularEdge[iedge][0]; - j=VerticesOfTriangularEdge[iedge][1]; - Real8 dij = detj-deti; - assert(i+j+k == 0 + 1 +2); - ba[j] = detj/dij; - ba[i] = -deti/dij; - ba[k] = 0; - ilast=NewItem(t,ba[0],ba[1],ba[2]); } - } // outside departure - - - - // recherche the intersection of [a,b] with Bh Mesh. - // we know a triangle ta contening the vertex a - // we have 2 case for intersection [a,b] with a edge [A,B] of Bh - // 1) the intersection point is in ]A,B[ - // 2) is A or B - // first version --- - for (;;) { - // t->Draw(); - if (iedge < 0) { - i0 =0;i1=1;i2=2; - dt[0] =bamg::det(a,b,(*t)[0]); - dt[1] =bamg::det(a,b,(*t)[1]); - dt[2] =bamg::det(a,b,(*t)[2]);} - else { - i2 = iedge; - i0 = NextEdge[i2]; - i1 = NextEdge[i0]; - dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because - dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side - dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);} - - // so we have just to see the transition from - to + of the det0..2 on edge of t - // because we are going from a to b - if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) && - ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0)) - ocut =i0; - else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) && - (dt[j=VerticesOfTriangularEdge[i1][1]] > 0)) - ocut =i1; - else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) && - (dt[j=VerticesOfTriangularEdge[i2][1]] > 0)) - ocut =i2; - else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) && - ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0)) - ocut =i0; - else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) && - (dt[j=VerticesOfTriangularEdge[i1][1]] > 0)) - ocut =i1; - else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) && - (dt[j=VerticesOfTriangularEdge[i2][1]] > 0)) - ocut =i2; - else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) && - ( dt[j=VerticesOfTriangularEdge[i0][1]] == 0)) - ocut =i0; - else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) && - (dt[j=VerticesOfTriangularEdge[i1][1]] == 0)) - ocut =i1; - else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) && - (dt[j=VerticesOfTriangularEdge[i2][1]] == 0)) - ocut =i2; - else { // On a edge (2 zero) - k =0; - if (dt[0]) ocut=0,k++; - if (dt[1]) ocut=1,k++; - if (dt[2]) ocut=2,k++; - if(k == 1) { - if (dt[ocut] >0) // triangle upper AB - ocut = NextEdge[ocut]; - i= VerticesOfTriangularEdge[ocut][0]; - j= VerticesOfTriangularEdge[ocut][1]; - } - else { - cerr << " Bug Split Edge " << endl; - cerr << " dt[0]= " << dt[0] - << " dt[1]= " << dt[1] - << " dt[2]= "<< dt[2] << endl; - cerr << i0 << " " << i1 << " " << i2 << endl; - cerr << " A = " << A << " B= " << B << endl; - cerr << " Triangle t = " << *t << endl; - cerr << (*t)[0] << (*t)[1] << (*t)[0] << endl; - cerr << " nbt = " << nbt << endl; - MeshError(100);}} - - k = OppositeVertex[ocut]; - - Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b); - - - if (detbij >= 0) { //we find the triangle contening b - dt[0]=bamg::det((*t)[1],(*t)[2],b); - dt[1]=bamg::det((*t)[2],(*t)[0],b); - dt[2]=bamg::det((*t)[0],(*t)[1],b); - Real8 dd = t->det; - NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd); - return ;} - else { // next triangle by adjacent by edge ocut - deti = dt[i]; - detj = dt[j]; - Real4 dij = detj-deti; - ba[i] = detj/dij; - ba[j] = -deti/dij; - ba[3-i-j ] = 0; - ilast=NewItem(t, ba[0],ba[1],ba[2]); - - TriangleAdjacent ta =t->Adj(ocut); - t = ta; - iedge= ta; - if (t->det <= 0) { - double ba,bb; - TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb); - NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1))); - // cout << " return " << ba << " " << bb << endl; - // ajoute le 03 frev 1997 par F. hecht - return; - } - }// we go outside of omega - } // for(;;) - - - } // routine SplitEdge - - - int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2) { - register int n; - R2 x(0,0); - if ( d0) x = (*tt)[0].r * d0; - if ( d1) x = x + (*tt)[1].r * d1; - if ( d2) x = x + (*tt)[2].r * d2; - // newer add same point - if(!Size || Norme2_2(lIntTria[Size-1].x-x)) { - if (Size==MaxSize) ReShape(); - lIntTria[Size].t=tt; - lIntTria[Size].bary[0]=d0; - lIntTria[Size].bary[1]=d1; - lIntTria[Size].bary[2]=d2; - lIntTria[Size].x = x; - Metric m0,m1,m2; - register Vertex * v; - if ((v=(*tt)(0))) m0 = v->m; - if ((v=(*tt)(1))) m1 = v->m; - if ((v=(*tt)(2))) m2 = v->m; - lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2); - n=Size++;} - else n=Size-1; - return n; - } - int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) { - register int n; - if(!Size || Norme2_2(lIntTria[Size-1].x-A)) { - if (Size==MaxSize) ReShape(); - lIntTria[Size].t=0; - lIntTria[Size].x=A; - lIntTria[Size].m=mm; - n=Size++; - } - else n=Size-1; - return n; - } - - Real8 ListofIntersectionTriangles::Length() - { - // cout << " n= " << Size << ":" ; - assert(Size>0); - // computation of the length - R2 C; - Metric Mx,My; - int ii,jj; - R2 x,y,xy; - - SegInterpolation *SegI=lSegsI; - SegI=lSegsI; - lSegsI[NbSeg].last=Size;// improvement - - int EndSeg=Size; - - y = lIntTria[0].x; - Real8 sxy, s = 0; - lIntTria[0].s =0; - SegI->lBegin=s; - - for (jj=0,ii=1;iilEnd=s, - SegI++, - EndSeg=SegI->last, - SegI->lBegin=s; - - // cout << ii << " " << jj << x<< y <lEnd=s; - - // cout << " len= " << s << endl; - return s; - } - - Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4 nbvx) - { - long int verbosity=0; - - - const Int4 nbvold = nbv; - Real8 s = Length(); - if (s < 1.5 ) return 0; - ////////////////////// - int ii = 1 ; - R2 y,x; - Metric My,Mx ; - Real8 sx =0,sy; - int nbi = Max(2,(int) (s+0.5)); - Real8 sint = s/nbi; - Real8 si = sint; - - int EndSeg=Size; - SegInterpolation *SegI=0; - if (NbSeg) - SegI=lSegsI,EndSeg=SegI->last; - - for (int k=1;klast; - - int ii1=ii-1; - x =lIntTria[ii1].x; - sx =lIntTria[ii1].s; - Metric Mx=lIntTria[ii1].m; - y =lIntTria[ii].x; - sy =lIntTria[ii].s; - Metric My=lIntTria[ii].m; - Real8 lxy = sy-sx; - Real8 cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy); - - R2 C; - Real8 cx = 1-cy; - C = SegI ? SegI->F(si): x * cx + y *cy; - - si += sint; - if ( nbvNumber(a) << ", " << CurrentTh->Number(b) << " ] and [ "; - cerr << CurrentTh->Number(aa) << " " << CurrentTh->Number(bb) << " ] " << endl; - } - MeshError(991); - */ } } Index: /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp (revision 2843) +++ /issm/trunk/src/c/Bamgx/objects/ListofIntersectionTriangles.cpp (revision 2843) @@ -0,0 +1,379 @@ +#include +#include +#include +#include + +#include "../Mesh2.h" +#include "../QuadTree.h" +#include "../SetOfE4.h" + +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../toolkits/toolkits.h" + +namespace bamg { + + /*Methods*/ + /*FUNCTION ListofIntersectionTriangles::SplitEdge{{{1*/ + void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh, const R2 &A,const R2 &B,int nbegin) { + Triangle *tbegin, *t; + + long int verbosity=2; + Icoor2 deta[3], deti,detj; + Real8 ba[3]; + int nbt =0,ifirst=-1,ilast; + int i0,i1,i2; + int ocut,i,j,k=-1; + // int OnAVertices =0; + Icoor2 dt[3]; + I2 a = Bh.toI2(A) ,b= Bh.toI2(B);// compute the Icoor a,b + I2 vi,vj; + int iedge =-1;// not a edge + + if(nbegin) {// optimisation + // we suppose knowing the starting triangle + t=tbegin=lIntTria[ilast=(Size-1)].t; + if (tbegin->det>=0) + ifirst = ilast;} + else {// not optimisation + init(); + t=tbegin = Bh.FindTriangleContening(a,deta); + if( t->det>=0) + ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det); + else + {// find the nearest boundary edge of the vertex A + // find a edge or such normal projection a the edge IJ is on the edge + // <=> IJ.IA >=0 && IJ.AJ >=0 + ilast=ifirst; + double ba,bb; + TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); + Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); + NewItem(A,Metric(ba,v0,bb,v1)); + t=edge; + // test if the point b is in the same side + if (det(v0.i,v1.i,b)>=0) { + //cout << " All the edge " << A << B << endl; + TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb); + Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1); + NewItem(A,Metric(ba,v0,bb,v1)); + return; + } + } // find the nearest boundary edge of the vertex A + } // end not optimisation + if (t->det<0) { // outside departure + while (t->det <0) { // intersection boundary edge and a,b, + k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; + assert(k>=0); + ocut = OppositeEdge[k]; + i=VerticesOfTriangularEdge[ocut][0]; + j=VerticesOfTriangularEdge[ocut][1]; + vi=(*t)[i]; + vj=(*t)[j]; + deti = bamg::det(a,b,vi); + detj = bamg::det(a,b,vj); + if (deti>0) // go to i direction on gamma + ocut = PreviousEdge[ocut]; + else if (detj<=0) // go to j direction on gamma + ocut = NextEdge[ocut]; + TriangleAdjacent tadj =t->Adj(ocut); + t = tadj; + iedge= tadj; + if (t == tbegin) { // + double ba,bb; + long int verbosity=2; + if (verbosity>7) + cout << " SplitEdge: All the edge " << A << B << nbegin << det(vi,vj,b) + << " deti= " << deti << " detj=" <det <0) + // theoriticaly we have: deti =<0 and detj>0 + + // computation of barycentric coor + // test if the point b is on size on t + // we revert vi,vj because vi,vj is def in Adj triangle + if ( det(vi,vj,b)>=0) { + if (verbosity>7) + cout << " SplitEdge: all AB outside " << A << B << endl; + t=tbegin; + Real8 ba,bb; + TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb); + NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1))); + return; + } + else + { + k = OppositeVertex[iedge]; + i=VerticesOfTriangularEdge[iedge][0]; + j=VerticesOfTriangularEdge[iedge][1]; + Real8 dij = detj-deti; + assert(i+j+k == 0 + 1 +2); + ba[j] = detj/dij; + ba[i] = -deti/dij; + ba[k] = 0; + ilast=NewItem(t,ba[0],ba[1],ba[2]); } + } // outside departure + + + + // recherche the intersection of [a,b] with Bh Mesh. + // we know a triangle ta contening the vertex a + // we have 2 case for intersection [a,b] with a edge [A,B] of Bh + // 1) the intersection point is in ]A,B[ + // 2) is A or B + // first version --- + for (;;) { + // t->Draw(); + if (iedge < 0) { + i0 =0;i1=1;i2=2; + dt[0] =bamg::det(a,b,(*t)[0]); + dt[1] =bamg::det(a,b,(*t)[1]); + dt[2] =bamg::det(a,b,(*t)[2]);} + else { + i2 = iedge; + i0 = NextEdge[i2]; + i1 = NextEdge[i0]; + dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because + dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side + dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);} + + // so we have just to see the transition from - to + of the det0..2 on edge of t + // because we are going from a to b + if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) && + ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0)) + ocut =i0; + else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i1][1]] > 0)) + ocut =i1; + else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i2][1]] > 0)) + ocut =i2; + else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) && + ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0)) + ocut =i0; + else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) && + (dt[j=VerticesOfTriangularEdge[i1][1]] > 0)) + ocut =i1; + else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) && + (dt[j=VerticesOfTriangularEdge[i2][1]] > 0)) + ocut =i2; + else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) && + ( dt[j=VerticesOfTriangularEdge[i0][1]] == 0)) + ocut =i0; + else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i1][1]] == 0)) + ocut =i1; + else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) && + (dt[j=VerticesOfTriangularEdge[i2][1]] == 0)) + ocut =i2; + else { // On a edge (2 zero) + k =0; + if (dt[0]) ocut=0,k++; + if (dt[1]) ocut=1,k++; + if (dt[2]) ocut=2,k++; + if(k == 1) { + if (dt[ocut] >0) // triangle upper AB + ocut = NextEdge[ocut]; + i= VerticesOfTriangularEdge[ocut][0]; + j= VerticesOfTriangularEdge[ocut][1]; + } + else { + cerr << " Bug Split Edge " << endl; + cerr << " dt[0]= " << dt[0] + << " dt[1]= " << dt[1] + << " dt[2]= "<< dt[2] << endl; + cerr << i0 << " " << i1 << " " << i2 << endl; + cerr << " A = " << A << " B= " << B << endl; + cerr << " Triangle t = " << *t << endl; + cerr << (*t)[0] << (*t)[1] << (*t)[0] << endl; + cerr << " nbt = " << nbt << endl; + MeshError(100);}} + + k = OppositeVertex[ocut]; + + Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b); + + + if (detbij >= 0) { //we find the triangle contening b + dt[0]=bamg::det((*t)[1],(*t)[2],b); + dt[1]=bamg::det((*t)[2],(*t)[0],b); + dt[2]=bamg::det((*t)[0],(*t)[1],b); + Real8 dd = t->det; + NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd); + return ;} + else { // next triangle by adjacent by edge ocut + deti = dt[i]; + detj = dt[j]; + Real4 dij = detj-deti; + ba[i] = detj/dij; + ba[j] = -deti/dij; + ba[3-i-j ] = 0; + ilast=NewItem(t, ba[0],ba[1],ba[2]); + + TriangleAdjacent ta =t->Adj(ocut); + t = ta; + iedge= ta; + if (t->det <= 0) { + double ba,bb; + TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb); + NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1))); + // cout << " return " << ba << " " << bb << endl; + // ajoute le 03 frev 1997 par F. hecht + return; + } + }// we go outside of omega + } // for(;;) + } + /*}}}1*/ + /*FUNCTION ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 disp('ok0');,Real8 disp('ok1');,Real8 d2) {{{1*/ + int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2) { + register int n; + R2 x(0,0); + if ( d0) x = (*tt)[0].r * d0; + if ( d1) x = x + (*tt)[1].r * d1; + if ( d2) x = x + (*tt)[2].r * d2; + // newer add same point + if(!Size || Norme2_2(lIntTria[Size-1].x-x)) { + if (Size==MaxSize) ReShape(); + lIntTria[Size].t=tt; + lIntTria[Size].bary[0]=d0; + lIntTria[Size].bary[1]=d1; + lIntTria[Size].bary[2]=d2; + lIntTria[Size].x = x; + Metric m0,m1,m2; + register Vertex * v; + if ((v=(*tt)(0))) m0 = v->m; + if ((v=(*tt)(1))) m1 = v->m; + if ((v=(*tt)(2))) m2 = v->m; + lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2); + n=Size++;} + else n=Size-1; + return n; + } + /*}}}1*/ + /*FUNCTION ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm){{{1*/ + int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) { + register int n; + if(!Size || Norme2_2(lIntTria[Size-1].x-A)) { + if (Size==MaxSize) ReShape(); + lIntTria[Size].t=0; + lIntTria[Size].x=A; + lIntTria[Size].m=mm; + n=Size++; + } + else n=Size-1; + return n; + } + /*}}}1*/ + /*FUNCTION ListofIntersectionTriangles::Length{{{1*/ + Real8 ListofIntersectionTriangles::Length(){ + // cout << " n= " << Size << ":" ; + assert(Size>0); + // computation of the length + R2 C; + Metric Mx,My; + int ii,jj; + R2 x,y,xy; + + SegInterpolation *SegI=lSegsI; + SegI=lSegsI; + lSegsI[NbSeg].last=Size;// improvement + + int EndSeg=Size; + + y = lIntTria[0].x; + Real8 sxy, s = 0; + lIntTria[0].s =0; + SegI->lBegin=s; + + for (jj=0,ii=1;iilEnd=s, + SegI++, + EndSeg=SegI->last, + SegI->lBegin=s; + + // cout << ii << " " << jj << x<< y <lEnd=s; + + // cout << " len= " << s << endl; + return s; + } + /*}}}1*/ + /*FUNCTION ListofIntersectionTriangles::NewPoints{{{1*/ + Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4 nbvx){ + long int verbosity=0; + + + const Int4 nbvold = nbv; + Real8 s = Length(); + if (s < 1.5 ) return 0; + ////////////////////// + int ii = 1 ; + R2 y,x; + Metric My,Mx ; + Real8 sx =0,sy; + int nbi = Max(2,(int) (s+0.5)); + Real8 sint = s/nbi; + Real8 si = sint; + + int EndSeg=Size; + SegInterpolation *SegI=0; + if (NbSeg) + SegI=lSegsI,EndSeg=SegI->last; + + for (int k=1;klast; + + int ii1=ii-1; + x =lIntTria[ii1].x; + sx =lIntTria[ii1].s; + Metric Mx=lIntTria[ii1].m; + y =lIntTria[ii].x; + sy =lIntTria[ii].s; + Metric My=lIntTria[ii].m; + Real8 lxy = sy-sx; + Real8 cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy); + + R2 C; + Real8 cx = 1-cy; + C = SegI ? SegI->F(si): x * cx + y *cy; + + si += sint; + if ( nbv s1 /___as2___\s2 ! + // \ a1|a2 / \ as1 / + // \ | / \ t1 / + // \ | / as2 \ a1/ + // \ | / \ / + // sa sa + // ------------------------------------------------------------- + int as1 = NextEdge[a1]; + int as2 = NextEdge[a2]; + int ap1 = PreviousEdge[a1]; + int ap2 = PreviousEdge[a2]; + (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb + (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa + // mise a jour des 2 adjacences externes + TriangleAdjacent taas1 = t1->Adj(as1), + taas2 = t2->Adj(as2), + tas1(t1,as1), tas2(t2,as2), + ta1(t1,a1),ta2(t2,a2); + // externe haut gauche + taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap()); + // externe bas droite + taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap()); + // remove the Mark UnMarkSwap + t1->SetUnMarkUnSwap(ap1); + t2->SetUnMarkUnSwap(ap2); + // interne + tas1.SetAdj2(tas2); + + t1->det = det1; + t2->det = det2; + + t1->SetTriangleContainingTheVertex(); + t2->SetTriangleContainingTheVertex(); + } // end swap + /*}}}1*/ + /*FUNCTION FindTriangle{{{1*/ + Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside){ + CurrentTh=&Th; + assert(&Th); + I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); + Icoor2 dete[3]; + Triangle & tb = *Th.FindTriangleContening(I,dete); + + if (tb.link) + { // internal point in a true triangles + a[0]= (Real8) dete[0]/ tb.det; + a[1]= (Real8) dete[1] / tb.det; + a[2] = (Real8) dete[2] / tb.det; + inside = 1; + return Th.Number(tb); + } + else + { + inside = 0; + double aa,bb; + TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb); + int k = ta; + Triangle * tc = ta; + if (!tc->link) + { ta = ta.Adj(); + tc=ta; + k = ta; + Exchange(aa,bb); + assert(tc->link); + } + a[VerticesOfTriangularEdge[k][0]] = aa; + a[VerticesOfTriangularEdge[k][1]] = bb; + a[OppositeVertex[k]] = 1- aa -bb; + return Th.Number(tc); + } + } + /*}}}1*/ + /*FUNCTION CloseBoundaryEdge{{{1*/ + TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) { + int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; + int dir=0; + assert(k>=0); + int kkk=0; + Icoor2 IJ_IA,IJ_AJ; + TriangleAdjacent edge(t,OppositeEdge[k]); + for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) + { + + assert(kkk++<1000); + Vertex &vI = *edge.EdgeVertex(0); + Vertex &vJ = *edge.EdgeVertex(1); + I2 I=vI, J=vJ, IJ= J-I; + IJ_IA = (IJ ,(A-I)); + // cout << A << vI.i << vJ.i << edge << " " << IJ_IA << " dir " << dir <0) {a=1;b=0;return edge;}// change of signe => I + else {dir=-1; + continue;}};// go in direction i + IJ_AJ = (IJ ,(J-A)); + if (IJ_AJ<0) { + if(dir<0) {a=0;b=1;return edge;} + else {dir = 1; + continue;}}// go in direction j + double IJ2 = IJ_IA + IJ_AJ; + assert(IJ2); + a= IJ_AJ/IJ2; + b= IJ_IA/IJ2; + // cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl; + return edge; + } + } + /*}}}1*/ + /*FUNCTION CloseBoundaryEdgeV2{{{1*/ + TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) { + // walk around the vertex + // version 2 for remove the probleme if we fill the hole + //int bug=1; + // Triangle *torigine = t; + // restart: + // int dir=0; + assert(t->link == 0); + // to have a starting edges + // try the 3 edge bourna-- in case of internal hole + // and choice the best + // + // the probleme is in case of the fine and long internal hole + // for exemple neart the training edge of a wing + Vertex * s=0,*s1=0, *s0=0; + Icoor2 imax = MaxICoor22; + Icoor2 l0 = imax,l1 = imax; + double dd2 = imax;// infinity + TriangleAdjacent er; + int cas=-2; + for (int j=0;j<3;j++) + { + TriangleAdjacent ta=t->FindBoundaryEdge(j); + if (! (Triangle *) ta) continue; + s0 = ta.EdgeVertex(0); + s1 = ta.EdgeVertex(1); + I2 A = * s0; + I2 B = *ta.EdgeVertex(1); + I2 AB = B-A,AC=C-A,BC=B-C; + Icoor2 ACAC = (AC,AC), BCBC = (BC,BC); + Icoor2 AB2 = Norme2_2(AB); // ||AB||^2 + Icoor2 ABAC = (AB,AC); // AB.AC| + + double d2; + if ( ABAC < 0 ) // DIST A + { + if ( (d2=(double) ACAC) < dd2) + { + // cout << " A " << d2 << " " << dd2; + er = ta; + l0 = ACAC; + l1 = BCBC; + cas = 0; + s = s0; + } + } + else if (ABAC > AB2) // DIST B + { + if ( (d2=(double) BCBC) < dd2) + { + // cout << " B " << d2 << " " << dd2; + dd2 = d2; + er = Adj(ta); // other direction + l0 = BCBC; + l1 = ACAC; + cas = 1; + s = s1; + } + } + else // DIST AB + { + + double det_2 = (double) Det(AB,AC); + det_2 *= det_2; // square of area*2 of triangle ABC + d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC + + if (d2 < dd2) + { + dd2 = d2; + er = ta; + l0 = (AC,AC); + l1 = (BC,BC); + s = 0; + cas = -1; + // cout << " ABAC " << ABAC << " ABAC " << ABAC + // << " AB2 " << AB2 << endl; + b = ((double) ABAC/(double) AB2); + a = 1 - b; + } + } + } + assert(cas !=-2); + // l1 = ||C s1|| , l0 = ||C s0|| + // where s0,s1 are the vertex of the edge er + + if ( s) + { + t=er; + TriangleAdjacent edge(er); + + int kkk=0; + int linkp = t->link == 0; + + Triangle * tt=t=edge=Adj(Previous(edge)); + do { // loop around vertex s + + assert(edge.EdgeVertex(0)==s && kkk++<10000); + + int link = tt->link == 0; + if ((link + linkp) == 1) + { // a boundary edge + Vertex * st = edge.EdgeVertex(1); + I2 I=*st; + Icoor2 ll = Norme2_2 (C-I); + if (ll < l1) { // the other vertex is neart + s1=st; + l1=ll; + er = edge; + if(ll