![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 2843
- Timestamp:
- 01/14/10 16:26:55 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 1 added
- 3 edited
-
Bamgx/Mesh2.cpp (modified) (5 diffs)
-
Bamgx/objects/ListofIntersectionTriangles.cpp (added)
-
Bamgx/objects/Triangles.cpp (modified) (1 diff)
-
Makefile.am (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Bamgx/Mesh2.cpp
r2814 r2843 1 // -*- Mode : c++ -*-2 //3 // SUMMARY :4 // USAGE :5 // ORG :6 // AUTHOR : Frederic Hecht7 // E-MAIL : hecht@ann.jussieu.fr8 //9 10 /*11 12 This file is part of Freefem++13 14 Freefem++ is free software; you can redistribute it and/or modify15 it under the terms of the GNU Lesser General Public License as published by16 the Free Software Foundation; either version 2.1 of the License, or17 (at your option) any later version.18 19 Freefem++ is distributed in the hope that it will be useful,20 but WITHOUT ANY WARRANTY; without even the implied warranty of21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the22 GNU Lesser General Public License for more details.23 24 You should have received a copy of the GNU Lesser General Public License25 along with Freefem++; if not, write to the Free Software26 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA27 */28 extern bool withrgraphique;29 30 1 #include "../shared/shared.h" 31 2 #include <stdio.h> … … 35 6 #include <iostream> 36 7 37 using namespace std;38 39 8 #include "Mesh2.h" 40 9 #include "QuadTree.h" 41 10 #include "SetOfE4.h" 42 11 12 using namespace std; 43 13 namespace bamg { 44 14 45 15 int Triangles::counter = 0; 46 47 Triangles * CurrentTh =0; 48 16 Triangles* CurrentTh=0; 49 17 int hinterpole=1; 50 51 18 int ForDebugging = 0; 52 19 const Direction NoDirOfSearch = Direction(); 53 #ifndef NDEBUG 54 inline void MyAssert(int i,char*ex,char * file,long line) 55 { 56 if( i) { 57 cerr << "Error Assert:" << ex << " in " << file << " line: " << line << endl; 58 #ifdef NOTFREEFEM 59 exit(1); 60 #else 61 throw(ErrorExec("exit",1000)); 62 #endif 63 } 64 } 65 #endif 66 67 Int4 AGoodNumberPrimeWith(Int4 n) 68 { 20 21 Int4 AGoodNumberPrimeWith(Int4 n){ 69 22 const Int4 BigPrimeNumber[] ={ 567890359L, 70 23 567890431L, 567890437L, 567890461L, 567890471L, … … 83 36 } 84 37 85 class Triangles;86 38 void MeshError(int Err,Triangles *Th){ 87 39 cerr << " Fatal error in the meshgenerator " << Err << endl ; 88 #ifdef NOTFREEFEM89 40 exit(1); 90 #else91 throw(ErrorMesh("Bamg",Err,Th));92 #endif93 41 } 94 42 95 ostream& operator <<(ostream& f, const Triangle & ta) 96 { 43 ostream& operator <<(ostream& f, const Triangle & ta) { 97 44 if(CurrentTh) 98 45 f << "[" << CurrentTh->Number(ta) << "::" … … 115 62 return f;} 116 63 117 void swap(Triangle *t1,Int1 a1,118 Triangle *t2,Int1 a2,119 Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2)120 { // swap121 // --------------------------------------------------------------122 // Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles123 //124 // sb sb125 // / | \ / \ !126 // as1/ | \ /a2 \ !127 // / | \ / t2 \ !128 // s1 /t1 | t2 \s2 --> s1 /___as2___\s2 !129 // \ a1|a2 / \ as1 /130 // \ | / \ t1 /131 // \ | / as2 \ a1/132 // \ | / \ /133 // sa sa134 // -------------------------------------------------------------135 int as1 = NextEdge[a1];136 int as2 = NextEdge[a2];137 int ap1 = PreviousEdge[a1];138 int ap2 = PreviousEdge[a2];139 (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb140 (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa141 // mise a jour des 2 adjacences externes142 TriangleAdjacent taas1 = t1->Adj(as1),143 taas2 = t2->Adj(as2),144 tas1(t1,as1), tas2(t2,as2),145 ta1(t1,a1),ta2(t2,a2);146 // externe haut gauche147 taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap());148 // externe bas droite149 taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap());150 // remove the Mark UnMarkSwap151 t1->SetUnMarkUnSwap(ap1);152 t2->SetUnMarkUnSwap(ap2);153 // interne154 tas1.SetAdj2(tas2);155 156 t1->det = det1;157 t2->det = det2;158 159 t1->SetTriangleContainingTheVertex();160 t2->SetTriangleContainingTheVertex();161 } // end swap162 163 164 165 166 167 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside)168 {169 CurrentTh=&Th;170 assert(&Th);171 I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y)));172 Icoor2 dete[3];173 Triangle & tb = *Th.FindTriangleContening(I,dete);174 175 if (tb.link)176 { // internal point in a true triangles177 a[0]= (Real8) dete[0]/ tb.det;178 a[1]= (Real8) dete[1] / tb.det;179 a[2] = (Real8) dete[2] / tb.det;180 inside = 1;181 return Th.Number(tb);182 }183 else184 {185 inside = 0;186 double aa,bb;187 TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb);188 int k = ta;189 Triangle * tc = ta;190 if (!tc->link)191 { ta = ta.Adj();192 tc=ta;193 k = ta;194 Exchange(aa,bb);195 assert(tc->link);196 }197 a[VerticesOfTriangularEdge[k][0]] = aa;198 a[VerticesOfTriangularEdge[k][1]] = bb;199 a[OppositeVertex[k]] = 1- aa -bb;200 return Th.Number(tc);201 }202 }203 204 205 TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) {206 //207 // cout << " - ";208 int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;209 int dir=0;210 assert(k>=0);211 int kkk=0;212 Icoor2 IJ_IA,IJ_AJ;213 TriangleAdjacent edge(t,OppositeEdge[k]);214 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge))))215 {216 217 assert(kkk++<1000);218 Vertex &vI = *edge.EdgeVertex(0);219 Vertex &vJ = *edge.EdgeVertex(1);220 I2 I=vI, J=vJ, IJ= J-I;221 IJ_IA = (IJ ,(A-I));222 // cout << A << vI.i << vJ.i << edge << " " << IJ_IA << " dir " << dir <<endl;223 if (IJ_IA<0) {224 if (dir>0) {a=1;b=0;return edge;}// change of signe => I225 else {dir=-1;226 continue;}};// go in direction i227 IJ_AJ = (IJ ,(J-A));228 if (IJ_AJ<0) {229 if(dir<0) {a=0;b=1;return edge;}230 else {dir = 1;231 continue;}}// go in direction j232 double IJ2 = IJ_IA + IJ_AJ;233 assert(IJ2);234 a= IJ_AJ/IJ2;235 b= IJ_IA/IJ2;236 // cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl;237 return edge;238 }239 }240 241 TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b)242 {243 // walk around the vertex244 // version 2 for remove the probleme if we fill the hole245 //int bug=1;246 // Triangle *torigine = t;247 // restart:248 // int dir=0;249 assert(t->link == 0);250 // to have a starting edges251 // try the 3 edge bourna-- in case of internal hole252 // and choice the best253 //254 //255 // the probleme is in case of the fine and long internal hole256 // for exemple neart the training edge of a wing257 //258 Vertex * s=0,*s1=0, *s0=0;259 Icoor2 imax = MaxICoor22;260 Icoor2 l0 = imax,l1 = imax;261 double dd2 = imax;// infinity262 TriangleAdjacent er;263 int cas=-2;264 for (int j=0;j<3;j++)265 {266 TriangleAdjacent ta=t->FindBoundaryEdge(j);267 if (! (Triangle *) ta) continue;268 s0 = ta.EdgeVertex(0);269 s1 = ta.EdgeVertex(1);270 I2 A = * s0;271 I2 B = *ta.EdgeVertex(1);272 I2 AB = B-A,AC=C-A,BC=B-C;273 Icoor2 ACAC = (AC,AC), BCBC = (BC,BC);274 Icoor2 AB2 = Norme2_2(AB); // ||AB||^2275 Icoor2 ABAC = (AB,AC); // AB.AC|276 277 double d2;278 if ( ABAC < 0 ) // DIST A279 {280 if ( (d2=(double) ACAC) < dd2)281 {282 // cout << " A " << d2 << " " << dd2;283 er = ta;284 l0 = ACAC;285 l1 = BCBC;286 cas = 0;287 s = s0;288 }289 }290 else if (ABAC > AB2) // DIST B291 {292 if ( (d2=(double) BCBC) < dd2)293 {294 // cout << " B " << d2 << " " << dd2;295 dd2 = d2;296 er = Adj(ta); // other direction297 l0 = BCBC;298 l1 = ACAC;299 cas = 1;300 s = s1;301 }302 }303 else // DIST AB304 {305 306 double det_2 = (double) Det(AB,AC);307 det_2 *= det_2; // square of area*2 of triangle ABC308 d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC309 // cout << " AB " << d2 << " " << dd2310 // << " " << CurrentTh->Number(ta.EdgeVertex(0))311 // << " " << CurrentTh->Number(ta.EdgeVertex(1)) << " " ;312 313 if (d2 < dd2)314 {315 dd2 = d2;316 er = ta;317 l0 = (AC,AC);318 l1 = (BC,BC);319 s = 0;320 cas = -1;321 // cout << " ABAC " << ABAC << " ABAC " << ABAC322 // << " AB2 " << AB2 << endl;323 b = ((double) ABAC/(double) AB2);324 a = 1 - b;325 }326 }327 }328 assert(cas !=-2);329 // l1 = ||C s1|| , l0 = ||C s0||330 // where s0,s1 are the vertex of the edge er331 332 if ( s)333 {334 t=er;335 TriangleAdjacent edge(er);336 337 int kkk=0;338 int linkp = t->link == 0;339 340 Triangle * tt=t=edge=Adj(Previous(edge));341 // cout << CurrentTh->Number(t) << " " << linkp << endl;342 do { // loop around vertex s343 344 assert(edge.EdgeVertex(0)==s && kkk++<10000);345 346 int link = tt->link == 0;347 // cout << CurrentTh->Number(tt) << " " << link << " " << CurrentTh->Number(s)348 // << " " << CurrentTh->Number(er.EdgeVertex(0))349 // << " " << CurrentTh->Number(er.EdgeVertex(1))350 // << " " << CurrentTh->Number(edge.EdgeVertex(0))351 // << " " << CurrentTh->Number(edge.EdgeVertex(1))352 // << endl;353 if ((link + linkp) == 1)354 { // a boundary edge355 Vertex * st = edge.EdgeVertex(1);356 I2 I=*st;357 Icoor2 ll = Norme2_2 (C-I);358 if (ll < l1) { // the other vertex is neart359 s1=st;360 l1=ll;361 er = edge;362 if(ll<l0) { // change of direction --363 s1=s;364 l1=l0;365 s=st;366 l0=ll;367 t=tt;368 edge=Adj(edge);369 link=linkp;370 er = edge;371 }372 }373 }374 375 linkp=link;376 edge=Adj(Previous(edge));377 tt = edge;378 } while (t!=tt);379 380 assert((Triangle *) er);381 I2 A((I2)*er.EdgeVertex(0));382 I2 B((I2)*er.EdgeVertex(1));383 I2 AB=B-A,AC=C-A,CB=B-C;384 double aa = (double) (AB,AC);385 double bb = (double) (AB,CB);386 // cout << " " << aa << " " << bb387 // << " " << CurrentTh->Number(er.EdgeVertex(0))388 // << " " << CurrentTh->Number(er.EdgeVertex(1)) ;389 if (aa<0) a=1,b=0;390 else if(bb<0) a=0,b=1;391 else392 {393 a = bb/(aa+bb);394 b = aa/(aa+bb);395 }396 }397 398 // cout <<" return= " << CurrentTh->Number(er.EdgeVertex(0)) << " "399 // << CurrentTh->Number(er.EdgeVertex(1)) << " " << a400 // << " " << b <<" " << l0 << " " <<l1 <<endl;401 return er;402 }403 404 405 void ListofIntersectionTriangles::SplitEdge(const Triangles & Bh,406 const R2 &A,const R2 &B,int nbegin)407 { // SplitEdge408 Triangle *tbegin, *t;409 410 long int verbosity=2;411 Icoor2 deta[3], deti,detj;412 Real8 ba[3];413 int nbt =0,ifirst=-1,ilast;414 int i0,i1,i2;415 int ocut,i,j,k=-1;416 // int OnAVertices =0;417 Icoor2 dt[3];418 I2 a = Bh.toI2(A) ,b= Bh.toI2(B);// compute the Icoor a,b419 I2 vi,vj;420 int iedge =-1;// not a edge421 422 if(nbegin) {// optimisation423 // we suppose knowing the starting triangle424 t=tbegin=lIntTria[ilast=(Size-1)].t;425 if (tbegin->det>=0)426 ifirst = ilast;}427 else {// not optimisation428 init();429 t=tbegin = Bh.FindTriangleContening(a,deta);430 if( t->det>=0)431 ilast=NewItem(t,Real8(deta[0])/t->det,Real8(deta[1])/t->det,Real8(deta[2])/t->det);432 else433 {// find the nearest boundary edge of the vertex A434 // find a edge or such normal projection a the edge IJ is on the edge435 // <=> IJ.IA >=0 && IJ.AJ >=0436 ilast=ifirst;437 double ba,bb;438 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);439 Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);440 NewItem(A,Metric(ba,v0,bb,v1));441 t=edge;442 // test if the point b is in the same side443 if (det(v0.i,v1.i,b)>=0) {444 //cout << " All the edge " << A << B << endl;445 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);446 Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);447 NewItem(A,Metric(ba,v0,bb,v1));448 return;449 }450 } // find the nearest boundary edge of the vertex A451 } // end not optimisation452 if (t->det<0) { // outside departure453 while (t->det <0) { // intersection boundary edge and a,b,454 k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0;455 assert(k>=0);456 ocut = OppositeEdge[k];457 i=VerticesOfTriangularEdge[ocut][0];458 j=VerticesOfTriangularEdge[ocut][1];459 vi=(*t)[i];460 vj=(*t)[j];461 deti = bamg::det(a,b,vi);462 detj = bamg::det(a,b,vj);463 if (deti>0) // go to i direction on gamma464 ocut = PreviousEdge[ocut];465 else if (detj<=0) // go to j direction on gamma466 ocut = NextEdge[ocut];467 TriangleAdjacent tadj =t->Adj(ocut);468 t = tadj;469 iedge= tadj;470 if (t == tbegin) { //471 double ba,bb;472 long int verbosity=2;473 if (verbosity>7)474 cout << " SplitEdge: All the edge " << A << B << nbegin << det(vi,vj,b)475 << " deti= " << deti << " detj=" <<detj << endl;476 TriangleAdjacent edge=CloseBoundaryEdge(a,t,ba,bb);477 Vertex & v0 = *edge.EdgeVertex(0), & v1 = *edge.EdgeVertex(1);478 NewItem(A,Metric(ba,v0,bb,v1));479 return;480 }481 } // end while (t->det <0)482 // theoriticaly we have: deti =<0 and detj>0483 484 // computation of barycentric coor485 // test if the point b is on size on t486 // we revert vi,vj because vi,vj is def in Adj triangle487 if ( det(vi,vj,b)>=0) {488 if (verbosity>7)489 cout << " SplitEdge: all AB outside " << A << B << endl;490 t=tbegin;491 Real8 ba,bb;492 TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);493 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));494 return;495 }496 else497 {498 k = OppositeVertex[iedge];499 i=VerticesOfTriangularEdge[iedge][0];500 j=VerticesOfTriangularEdge[iedge][1];501 Real8 dij = detj-deti;502 assert(i+j+k == 0 + 1 +2);503 ba[j] = detj/dij;504 ba[i] = -deti/dij;505 ba[k] = 0;506 ilast=NewItem(t,ba[0],ba[1],ba[2]); }507 } // outside departure508 509 510 511 // recherche the intersection of [a,b] with Bh Mesh.512 // we know a triangle ta contening the vertex a513 // we have 2 case for intersection [a,b] with a edge [A,B] of Bh514 // 1) the intersection point is in ]A,B[515 // 2) is A or B516 // first version ---517 for (;;) {518 // t->Draw();519 if (iedge < 0) {520 i0 =0;i1=1;i2=2;521 dt[0] =bamg::det(a,b,(*t)[0]);522 dt[1] =bamg::det(a,b,(*t)[1]);523 dt[2] =bamg::det(a,b,(*t)[2]);}524 else {525 i2 = iedge;526 i0 = NextEdge[i2];527 i1 = NextEdge[i0];528 dt[VerticesOfTriangularEdge[iedge][0]] = detj;// we revert i,j because529 dt[VerticesOfTriangularEdge[iedge][1]] = deti;// we take the Triangle by the other side530 dt[iedge] = det(a,b,(*t)[OppositeVertex[iedge]]);}531 532 // so we have just to see the transition from - to + of the det0..2 on edge of t533 // because we are going from a to b534 if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&535 ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))536 ocut =i0;537 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&538 (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))539 ocut =i1;540 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&541 (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))542 ocut =i2;543 else if ((dt[i=VerticesOfTriangularEdge[i0][0]] == 0) &&544 ( dt[j=VerticesOfTriangularEdge[i0][1]] > 0))545 ocut =i0;546 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] == 0) &&547 (dt[j=VerticesOfTriangularEdge[i1][1]] > 0))548 ocut =i1;549 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] == 0) &&550 (dt[j=VerticesOfTriangularEdge[i2][1]] > 0))551 ocut =i2;552 else if ((dt[i=VerticesOfTriangularEdge[i0][0]] < 0) &&553 ( dt[j=VerticesOfTriangularEdge[i0][1]] == 0))554 ocut =i0;555 else if ((dt[i=VerticesOfTriangularEdge[i1][0]] < 0) &&556 (dt[j=VerticesOfTriangularEdge[i1][1]] == 0))557 ocut =i1;558 else if ((dt[i=VerticesOfTriangularEdge[i2][0]] < 0) &&559 (dt[j=VerticesOfTriangularEdge[i2][1]] == 0))560 ocut =i2;561 else { // On a edge (2 zero)562 k =0;563 if (dt[0]) ocut=0,k++;564 if (dt[1]) ocut=1,k++;565 if (dt[2]) ocut=2,k++;566 if(k == 1) {567 if (dt[ocut] >0) // triangle upper AB568 ocut = NextEdge[ocut];569 i= VerticesOfTriangularEdge[ocut][0];570 j= VerticesOfTriangularEdge[ocut][1];571 }572 else {573 cerr << " Bug Split Edge " << endl;574 cerr << " dt[0]= " << dt[0]575 << " dt[1]= " << dt[1]576 << " dt[2]= "<< dt[2] << endl;577 cerr << i0 << " " << i1 << " " << i2 << endl;578 cerr << " A = " << A << " B= " << B << endl;579 cerr << " Triangle t = " << *t << endl;580 cerr << (*t)[0] << (*t)[1] << (*t)[0] << endl;581 cerr << " nbt = " << nbt << endl;582 MeshError(100);}}583 584 k = OppositeVertex[ocut];585 586 Icoor2 detbij = bamg::det((*t)[i],(*t)[j],b);587 588 589 if (detbij >= 0) { //we find the triangle contening b590 dt[0]=bamg::det((*t)[1],(*t)[2],b);591 dt[1]=bamg::det((*t)[2],(*t)[0],b);592 dt[2]=bamg::det((*t)[0],(*t)[1],b);593 Real8 dd = t->det;594 NewItem(t,dt[0]/dd,dt[1]/dd,dt[2]/dd);595 return ;}596 else { // next triangle by adjacent by edge ocut597 deti = dt[i];598 detj = dt[j];599 Real4 dij = detj-deti;600 ba[i] = detj/dij;601 ba[j] = -deti/dij;602 ba[3-i-j ] = 0;603 ilast=NewItem(t, ba[0],ba[1],ba[2]);604 605 TriangleAdjacent ta =t->Adj(ocut);606 t = ta;607 iedge= ta;608 if (t->det <= 0) {609 double ba,bb;610 TriangleAdjacent edge=CloseBoundaryEdge(b,t,ba,bb);611 NewItem(B,Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1)));612 // cout << " return " << ba << " " << bb << endl;613 // ajoute le 03 frev 1997 par F. hecht614 return;615 }616 }// we go outside of omega617 } // for(;;)618 619 620 } // routine SplitEdge621 622 623 int ListofIntersectionTriangles::NewItem(Triangle * tt,Real8 d0,Real8 d1,Real8 d2) {624 register int n;625 R2 x(0,0);626 if ( d0) x = (*tt)[0].r * d0;627 if ( d1) x = x + (*tt)[1].r * d1;628 if ( d2) x = x + (*tt)[2].r * d2;629 // newer add same point630 if(!Size || Norme2_2(lIntTria[Size-1].x-x)) {631 if (Size==MaxSize) ReShape();632 lIntTria[Size].t=tt;633 lIntTria[Size].bary[0]=d0;634 lIntTria[Size].bary[1]=d1;635 lIntTria[Size].bary[2]=d2;636 lIntTria[Size].x = x;637 Metric m0,m1,m2;638 register Vertex * v;639 if ((v=(*tt)(0))) m0 = v->m;640 if ((v=(*tt)(1))) m1 = v->m;641 if ((v=(*tt)(2))) m2 = v->m;642 lIntTria[Size].m = Metric(lIntTria[Size].bary,m0,m1,m2);643 n=Size++;}644 else n=Size-1;645 return n;646 }647 int ListofIntersectionTriangles::NewItem(R2 A,const Metric & mm) {648 register int n;649 if(!Size || Norme2_2(lIntTria[Size-1].x-A)) {650 if (Size==MaxSize) ReShape();651 lIntTria[Size].t=0;652 lIntTria[Size].x=A;653 lIntTria[Size].m=mm;654 n=Size++;655 }656 else n=Size-1;657 return n;658 }659 660 Real8 ListofIntersectionTriangles::Length()661 {662 // cout << " n= " << Size << ":" ;663 assert(Size>0);664 // computation of the length665 R2 C;666 Metric Mx,My;667 int ii,jj;668 R2 x,y,xy;669 670 SegInterpolation *SegI=lSegsI;671 SegI=lSegsI;672 lSegsI[NbSeg].last=Size;// improvement673 674 int EndSeg=Size;675 676 y = lIntTria[0].x;677 Real8 sxy, s = 0;678 lIntTria[0].s =0;679 SegI->lBegin=s;680 681 for (jj=0,ii=1;ii<Size;jj=ii++)682 {683 // seg jj,ii684 x=y;685 y = lIntTria[ii].x;686 xy = y-x;687 Mx = lIntTria[ii].m;688 My = lIntTria[jj].m;689 // Real8 &sx= lIntTria[ii].sp; // previous seg690 // Real8 &sy= lIntTria[jj].sn; // next seg691 // sx = Mx(xy);692 // sy = My(xy);693 // sxy = (Mx(xy)+ My(xy))/2.0;694 sxy = LengthInterpole(Mx,My,xy);695 s += sxy;696 lIntTria[ii].s = s;697 if (ii == EndSeg)698 SegI->lEnd=s,699 SegI++,700 EndSeg=SegI->last,701 SegI->lBegin=s;702 703 // cout << ii << " " << jj << x<< y <<xy << s << lIntTria[ii].m ;704 }705 len = s;706 SegI->lEnd=s;707 708 // cout << " len= " << s << endl;709 return s;710 }711 712 Int4 ListofIntersectionTriangles::NewPoints(Vertex * vertices,Int4 & nbv,Int4 nbvx)713 {714 long int verbosity=0;715 716 717 const Int4 nbvold = nbv;718 Real8 s = Length();719 if (s < 1.5 ) return 0;720 //////////////////////721 int ii = 1 ;722 R2 y,x;723 Metric My,Mx ;724 Real8 sx =0,sy;725 int nbi = Max(2,(int) (s+0.5));726 Real8 sint = s/nbi;727 Real8 si = sint;728 729 int EndSeg=Size;730 SegInterpolation *SegI=0;731 if (NbSeg)732 SegI=lSegsI,EndSeg=SegI->last;733 734 for (int k=1;k<nbi;k++)735 {736 while ((ii < Size) && ( lIntTria[ii].s <= si ))737 if (ii++ == EndSeg)738 SegI++,EndSeg=SegI->last;739 740 int ii1=ii-1;741 x =lIntTria[ii1].x;742 sx =lIntTria[ii1].s;743 Metric Mx=lIntTria[ii1].m;744 y =lIntTria[ii].x;745 sy =lIntTria[ii].s;746 Metric My=lIntTria[ii].m;747 Real8 lxy = sy-sx;748 Real8 cy = abscisseInterpole(Mx,My,y-x,(si-sx)/lxy);749 750 R2 C;751 Real8 cx = 1-cy;752 C = SegI ? SegI->F(si): x * cx + y *cy;753 754 si += sint;755 if ( nbv<nbvx) {756 vertices[nbv].r = C;757 vertices[nbv++].m = Metric(cx,lIntTria[ii-1].m,cy,lIntTria[ii].m);758 if((verbosity/100%10)==2)759 cout << " -- Add point " << nbv-1 << " " << vertices[nbv-1] << " " << vertices[nbv-1].m << endl;760 761 }762 else return nbv-nbvold;763 }764 return nbv-nbvold;765 }766 64 767 65 int SwapForForcingEdge(Vertex * & pva ,Vertex * & pvb , … … 908 206 taret = tc; 909 207 return -2; // error boundary is crossing 910 /* cerr << "Fatal Error boundary is crossing ";911 if(CurrentTh)912 {913 cerr << " edge: [" << CurrentTh->Number(a) << ", " << CurrentTh->Number(b) << " ] and [ ";914 cerr << CurrentTh->Number(aa) << " " << CurrentTh->Number(bb) << " ] " << endl;915 }916 MeshError(991);917 */918 208 } 919 209 } -
issm/trunk/src/c/Bamgx/objects/Triangles.cpp
r2841 r2843 5762 5762 /*}}}1*/ 5763 5763 5764 /*Intermediary*/ 5765 /*FUNCTION swap{{{1*/ 5766 void swap(Triangle *t1,Int1 a1, Triangle *t2,Int1 a2, Vertex *s1,Vertex *s2,Icoor2 det1,Icoor2 det2){ 5767 // -------------------------------------------------------------- 5768 // Int1 a2=aa[a];// les 2 numero de l arete dans les 2 triangles 5769 // 5770 // sb sb 5771 // / | \ / \ ! 5772 // as1/ | \ /a2 \ ! 5773 // / | \ / t2 \ ! 5774 // s1 /t1 | t2 \s2 --> s1 /___as2___\s2 ! 5775 // \ a1|a2 / \ as1 / 5776 // \ | / \ t1 / 5777 // \ | / as2 \ a1/ 5778 // \ | / \ / 5779 // sa sa 5780 // ------------------------------------------------------------- 5781 int as1 = NextEdge[a1]; 5782 int as2 = NextEdge[a2]; 5783 int ap1 = PreviousEdge[a1]; 5784 int ap2 = PreviousEdge[a2]; 5785 (*t1)(VerticesOfTriangularEdge[a1][1]) = s2 ; // avant sb 5786 (*t2)(VerticesOfTriangularEdge[a2][1]) = s1 ; // avant sa 5787 // mise a jour des 2 adjacences externes 5788 TriangleAdjacent taas1 = t1->Adj(as1), 5789 taas2 = t2->Adj(as2), 5790 tas1(t1,as1), tas2(t2,as2), 5791 ta1(t1,a1),ta2(t2,a2); 5792 // externe haut gauche 5793 taas1.SetAdj2(ta2, taas1.GetAllFlag_UnSwap()); 5794 // externe bas droite 5795 taas2.SetAdj2(ta1, taas2.GetAllFlag_UnSwap()); 5796 // remove the Mark UnMarkSwap 5797 t1->SetUnMarkUnSwap(ap1); 5798 t2->SetUnMarkUnSwap(ap2); 5799 // interne 5800 tas1.SetAdj2(tas2); 5801 5802 t1->det = det1; 5803 t2->det = det2; 5804 5805 t1->SetTriangleContainingTheVertex(); 5806 t2->SetTriangleContainingTheVertex(); 5807 } // end swap 5808 /*}}}1*/ 5809 /*FUNCTION FindTriangle{{{1*/ 5810 Int4 FindTriangle(Triangles &Th, Real8 x, Real8 y, double* a,int & inside){ 5811 CurrentTh=&Th; 5812 assert(&Th); 5813 I2 I = Th.toI2(R2(Min(Max(Th.pmin.x,x),Th.pmax.x),Min(Max(Th.pmin.y,y),Th.pmax.y))); 5814 Icoor2 dete[3]; 5815 Triangle & tb = *Th.FindTriangleContening(I,dete); 5816 5817 if (tb.link) 5818 { // internal point in a true triangles 5819 a[0]= (Real8) dete[0]/ tb.det; 5820 a[1]= (Real8) dete[1] / tb.det; 5821 a[2] = (Real8) dete[2] / tb.det; 5822 inside = 1; 5823 return Th.Number(tb); 5824 } 5825 else 5826 { 5827 inside = 0; 5828 double aa,bb; 5829 TriangleAdjacent ta=CloseBoundaryEdgeV2(I,&tb,aa,bb); 5830 int k = ta; 5831 Triangle * tc = ta; 5832 if (!tc->link) 5833 { ta = ta.Adj(); 5834 tc=ta; 5835 k = ta; 5836 Exchange(aa,bb); 5837 assert(tc->link); 5838 } 5839 a[VerticesOfTriangularEdge[k][0]] = aa; 5840 a[VerticesOfTriangularEdge[k][1]] = bb; 5841 a[OppositeVertex[k]] = 1- aa -bb; 5842 return Th.Number(tc); 5843 } 5844 } 5845 /*}}}1*/ 5846 /*FUNCTION CloseBoundaryEdge{{{1*/ 5847 TriangleAdjacent CloseBoundaryEdge(I2 A,Triangle *t, double &a,double &b) { 5848 int k=(*t)(0) ? (( (*t)(1) ? ( (*t)(2) ? -1 : 2) : 1 )) : 0; 5849 int dir=0; 5850 assert(k>=0); 5851 int kkk=0; 5852 Icoor2 IJ_IA,IJ_AJ; 5853 TriangleAdjacent edge(t,OppositeEdge[k]); 5854 for (;;edge = dir >0 ? Next(Adj(Next(edge))) : Previous(Adj(Previous(edge)))) 5855 { 5856 5857 assert(kkk++<1000); 5858 Vertex &vI = *edge.EdgeVertex(0); 5859 Vertex &vJ = *edge.EdgeVertex(1); 5860 I2 I=vI, J=vJ, IJ= J-I; 5861 IJ_IA = (IJ ,(A-I)); 5862 // cout << A << vI.i << vJ.i << edge << " " << IJ_IA << " dir " << dir <<endl; 5863 if (IJ_IA<0) { 5864 if (dir>0) {a=1;b=0;return edge;}// change of signe => I 5865 else {dir=-1; 5866 continue;}};// go in direction i 5867 IJ_AJ = (IJ ,(J-A)); 5868 if (IJ_AJ<0) { 5869 if(dir<0) {a=0;b=1;return edge;} 5870 else {dir = 1; 5871 continue;}}// go in direction j 5872 double IJ2 = IJ_IA + IJ_AJ; 5873 assert(IJ2); 5874 a= IJ_AJ/IJ2; 5875 b= IJ_IA/IJ2; 5876 // cout<< "CloseBoundaryEdge a = " << a << " b= " << b << endl; 5877 return edge; 5878 } 5879 } 5880 /*}}}1*/ 5881 /*FUNCTION CloseBoundaryEdgeV2{{{1*/ 5882 TriangleAdjacent CloseBoundaryEdgeV2(I2 C,Triangle *t, double &a,double &b) { 5883 // walk around the vertex 5884 // version 2 for remove the probleme if we fill the hole 5885 //int bug=1; 5886 // Triangle *torigine = t; 5887 // restart: 5888 // int dir=0; 5889 assert(t->link == 0); 5890 // to have a starting edges 5891 // try the 3 edge bourna-- in case of internal hole 5892 // and choice the best 5893 // 5894 // the probleme is in case of the fine and long internal hole 5895 // for exemple neart the training edge of a wing 5896 Vertex * s=0,*s1=0, *s0=0; 5897 Icoor2 imax = MaxICoor22; 5898 Icoor2 l0 = imax,l1 = imax; 5899 double dd2 = imax;// infinity 5900 TriangleAdjacent er; 5901 int cas=-2; 5902 for (int j=0;j<3;j++) 5903 { 5904 TriangleAdjacent ta=t->FindBoundaryEdge(j); 5905 if (! (Triangle *) ta) continue; 5906 s0 = ta.EdgeVertex(0); 5907 s1 = ta.EdgeVertex(1); 5908 I2 A = * s0; 5909 I2 B = *ta.EdgeVertex(1); 5910 I2 AB = B-A,AC=C-A,BC=B-C; 5911 Icoor2 ACAC = (AC,AC), BCBC = (BC,BC); 5912 Icoor2 AB2 = Norme2_2(AB); // ||AB||^2 5913 Icoor2 ABAC = (AB,AC); // AB.AC| 5914 5915 double d2; 5916 if ( ABAC < 0 ) // DIST A 5917 { 5918 if ( (d2=(double) ACAC) < dd2) 5919 { 5920 // cout << " A " << d2 << " " << dd2; 5921 er = ta; 5922 l0 = ACAC; 5923 l1 = BCBC; 5924 cas = 0; 5925 s = s0; 5926 } 5927 } 5928 else if (ABAC > AB2) // DIST B 5929 { 5930 if ( (d2=(double) BCBC) < dd2) 5931 { 5932 // cout << " B " << d2 << " " << dd2; 5933 dd2 = d2; 5934 er = Adj(ta); // other direction 5935 l0 = BCBC; 5936 l1 = ACAC; 5937 cas = 1; 5938 s = s1; 5939 } 5940 } 5941 else // DIST AB 5942 { 5943 5944 double det_2 = (double) Det(AB,AC); 5945 det_2 *= det_2; // square of area*2 of triangle ABC 5946 d2 = det_2/ (double) AB2; // hauteur^2 in C of of triangle ABC 5947 5948 if (d2 < dd2) 5949 { 5950 dd2 = d2; 5951 er = ta; 5952 l0 = (AC,AC); 5953 l1 = (BC,BC); 5954 s = 0; 5955 cas = -1; 5956 // cout << " ABAC " << ABAC << " ABAC " << ABAC 5957 // << " AB2 " << AB2 << endl; 5958 b = ((double) ABAC/(double) AB2); 5959 a = 1 - b; 5960 } 5961 } 5962 } 5963 assert(cas !=-2); 5964 // l1 = ||C s1|| , l0 = ||C s0|| 5965 // where s0,s1 are the vertex of the edge er 5966 5967 if ( s) 5968 { 5969 t=er; 5970 TriangleAdjacent edge(er); 5971 5972 int kkk=0; 5973 int linkp = t->link == 0; 5974 5975 Triangle * tt=t=edge=Adj(Previous(edge)); 5976 do { // loop around vertex s 5977 5978 assert(edge.EdgeVertex(0)==s && kkk++<10000); 5979 5980 int link = tt->link == 0; 5981 if ((link + linkp) == 1) 5982 { // a boundary edge 5983 Vertex * st = edge.EdgeVertex(1); 5984 I2 I=*st; 5985 Icoor2 ll = Norme2_2 (C-I); 5986 if (ll < l1) { // the other vertex is neart 5987 s1=st; 5988 l1=ll; 5989 er = edge; 5990 if(ll<l0) { // change of direction -- 5991 s1=s; 5992 l1=l0; 5993 s=st; 5994 l0=ll; 5995 t=tt; 5996 edge=Adj(edge); 5997 link=linkp; 5998 er = edge; 5999 } 6000 } 6001 } 6002 6003 linkp=link; 6004 edge=Adj(Previous(edge)); 6005 tt = edge; 6006 } while (t!=tt); 6007 6008 assert((Triangle *) er); 6009 I2 A((I2)*er.EdgeVertex(0)); 6010 I2 B((I2)*er.EdgeVertex(1)); 6011 I2 AB=B-A,AC=C-A,CB=B-C; 6012 double aa = (double) (AB,AC); 6013 double bb = (double) (AB,CB); 6014 if (aa<0) a=1,b=0; 6015 else if(bb<0) a=0,b=1; 6016 else 6017 { 6018 a = bb/(aa+bb); 6019 b = aa/(aa+bb); 6020 } 6021 } 6022 return er; 6023 } 6024 /*}}}1*/ 6025 5764 6026 } -
issm/trunk/src/c/Makefile.am
r2842 r2843 326 326 ./Bamgx/objects/MetricAnIso.cpp \ 327 327 ./Bamgx/objects/GeometricalEdge.cpp \ 328 ./Bamgx/objects/ListofIntersectionTriangles.cpp \ 328 329 ./Bamgx/objects/Triangles.cpp \ 329 330 ./Bamgx/objects/Triangle.cpp \ … … 666 667 ./Bamgx/objects/MetricAnIso.cpp \ 667 668 ./Bamgx/objects/GeometricalEdge.cpp \ 669 ./Bamgx/objects/ListofIntersectionTriangles.cpp \ 668 670 ./Bamgx/objects/Triangles.cpp \ 669 671 ./Bamgx/objects/Triangle.cpp \
Note:
See TracChangeset
for help on using the changeset viewer.
